Why bvp4/5c produces anomalous solution when using parameter passing using anonymous functions ?

I've found the following strange behaviour of the bvp solvers:

taking the example in matlab "shockbvp" and modifying it so that instead of nested functions it uses anonymous functions for passing the parameter "e" in the example, both solvers bvp4/5c produce huge errors when "e" decreases to 10^(-3).

Maybe I've made some stupid mistake...

I'm under Win10, using Matlab 2017b.

See the modified shockbvp code file below :

%function shockbvp(solver)
%SHOCKBVP  The solution has a shock layer near x = 0
%   This is an example used in U. Ascher, R. Mattheij, and R. Russell,
%   Numerical Solution of Boundary Value Problems for Ordinary Differential
%   Equations, SIAM, Philadelphia, PA, 1995,  to illustrate the mesh
%   selection strategy of COLSYS.
%
%   For 0 < e << 1, the solution of
%
%       e*y'' + x*y' = -e*pi^2*cos(pi*x) - pi*x*sin(pi*x)
%
%   on the interval [-1,1] with boundary conditions y(-1) = -2 and y(1) = 0
%   has a rapid transition layer at x = 0.
%
%   This example illustrates how a numerically difficult problem (e = 1e-4)
%   can be solved successfully using continuation. For this problem,
%   analytical partial derivatives are easy to derive and the solver benefits
%   from using them.
%
%   By default, this example uses the BVP4C solver. Use syntax
%   SHOCKBVP('bvp5c') to solve this problem with the BVP5C solver.
%
%   See also BVP4C, BVP5C, BVPSET, BVPGET, BVPINIT, DEVAL, FUNCTION_HANDLE.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

%if nargin < 1
%   solver = 'bvp4c';
%end
%bvpsolver = fcnchk(solver);

% The differential equations written as a first order system and the
% boundary conditions are coded in shockODE and shockBC, respectively. Their
% partial derivatives are coded in shockJac and shockBCJac and passed to the
% solver via the options. The option 'Vectorized' instructs the solver that
% the differential equation function has been vectorized, i.e.
% shockODE([x1 x2 ...],[y1 y2 ...]) returns [shockODE(x1,y1) shockODE(x2,y2) ...].
% Such coding improves the solver performance.

e=0.1;
options = bvpset('FJacobian',@(x,y)shockJac(x,y,e),'BCJacobian',@(x,y)shockBCJac(x,y,e),'Vectorized','on','stats','on');
%options = bvpset('FJacobian',@(x,y)shockJac(x,y,e),'BCJacobian',@(x,y)shockBCJac(x,y,e),'Vectorized','on','stats','on','RelTol', 1e-10,'AbsTol',1e-11,'Nmax',15000);
%options = bvpset('FJacobian',@shockJac,'BCJacobian',@shockBCJac,'Vectorized','on');
% A guess for the initial mesh and the solution
sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

% Problem parameter e is shared with the nested functions.
% Solution for e = 1e-2, 1e-3, 1e-4 obtained using continuation.
%e = 0.1;
for i = 2:2
   e = e/10;
   sol = bvp5c(@(x,y)shockODE(x,y,e),@(x,y)shockBC(x,y,e),sol,options);
end

Warning: Unable to meet the tolerance without using more than 5000 mesh points.
The last mesh of 1117 points and the solution are available in the output argument.
The maximum error is 4.38728, while requested accuracy is 0.001.

The solution was obtained on a mesh of 1117 points.
The maximum error is  4.387e+00. 
There were 3997 calls to the ODE function. 
There were 123 calls to the BC function. 

%sol = bvp4c(@(x,y)shockODE(x,y,e),@(x,y)shockBC(x,y,e),sol,options);
% The final solution
figure;
plot(sol.x,sol.y(1,:));
axis([-1 1 -2.2 2.2]);
title(['There is a shock at x = 0 when \epsilon =' sprintf('%.e',e) '.']);
xlabel('x');
ylabel('solution y');

% -----------------------------------------------------------------------
% Nested functions -- e is shared with the outer function.
%

   function dydx = shockODE(x,y,e)
      % SHOCKODE  Evaluate the ODE function (vectorized)
      pix = pi*x;
      dydx = [                 y(2,:)
         -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix) ];
   end
% -----------------------------------------------------------------------

   function res = shockBC(ya,yb,e)
      % SHOCKBC  Evaluate the residual in the boundary conditions
      res = [ ya(1)+2
         yb(1)  ];
   end
% -----------------------------------------------------------------------

   function jac = shockJac(x,y,e)
      % SHOCKJAC  Evaluate the Jacobian of the ODE function
      %  x and y are required arguments.
      jac = [ 0   1
         0 -x/e ];
   end
% -----------------------------------------------------------------------

   function [dBCdya,dBCdyb] = shockBCJac(ya,yb,e)
      % SHOCKBCJAC  Evaluate the partial derivatives of the boundary conditions
      %  ya and yb are required arguments.
      dBCdya = [ 1 0
         0 0 ];
      
      dBCdyb = [ 0 0
         1 0 ];
   end
% -----------------------------------------------------------------------

%end  % shockbvp

Answer 1

[Torsten]

Works for me.

With each call to bvp5c, the options structure has to be modified because e changes within the loop.

e = 1;
sol = bvpinit(linspace(-1,1,1000),[1 0]);
hold on
for i = 1:6
   options = bvpset('FJacobian',@(x,y)shockJac(x,y,e),'BCJacobian',@(x,y)shockBCJac(x,y,e),'Vectorized','on','stats','on');
   sol = bvp5c(@(x,y)shockODE(x,y,e),@(x,y)shockBC(x,y,e),sol,options);
   e = e/10;
   plot(sol.x,sol.y(1,:));   
end

The solution was obtained on a mesh of 1000 points.
The maximum error is  1.240e-11. 
There were 2001 calls to the ODE function. 
There were 3 calls to the BC function. 
The solution was obtained on a mesh of 1000 points.
The maximum error is  3.250e-11. 
There were 2001 calls to the ODE function. 
There were 3 calls to the BC function. 
The solution was obtained on a mesh of 1000 points.
The maximum error is  3.387e-09. 
There were 2001 calls to the ODE function. 
There were 3 calls to the BC function. 
The solution was obtained on a mesh of 1000 points.
The maximum error is  1.030e-06. 
There were 2001 calls to the ODE function. 
There were 3 calls to the BC function. 
The solution was obtained on a mesh of 1000 points.
The maximum error is  3.094e-04. 
There were 2001 calls to the ODE function. 
There were 3 calls to the BC function. 
The solution was obtained on a mesh of 323 points.
The maximum error is  5.695e-04. 
There were 3635 calls to the ODE function. 
There were 9 calls to the BC function. 

hold off
axis([-1 1 -2.2 2.2]);
xlabel('x');
ylabel('solution y')
grid on

function dydx = shockODE(x,y,e)
  % SHOCKODE  Evaluate the ODE function (vectorized)
  pix = pi*x;
  dydx = [y(2,:)
         -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix) ];
end
% -----------------------------------------------------------------------
function res = shockBC(ya,yb,e)
  % SHOCKBC  Evaluate the residual in the boundary conditions
  res = [ ya(1)+2
          yb(1)  ];
end
% -----------------------------------------------------------------------
function jac = shockJac(x,y,e)
   % SHOCKJAC  Evaluate the Jacobian of the ODE function
   %  x and y are required arguments.
   jac = [ 0   1
           0 -x/e ];
end
% -----------------------------------------------------------------------
function [dBCdya,dBCdyb] = shockBCJac(ya,yb,e)
   % SHOCKBCJAC  Evaluate the partial derivatives of the boundary conditions
   %  ya and yb are required arguments.
   dBCdya = [ 1 0
              0 0 ];
      
   dBCdyb = [ 0 0
              1 0 ];
end

[Peter Lefebre]

Thanks Torsten for this enlightening answer!

It's always good to learn.

Best, Peter