gusucode.com > demos工具箱matlab源码程序 > demos/shockbvp.m

    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.

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:4
   e = e/10;
   sol = bvpsolver(@shockODE,@shockBC,sol,options);
end

% 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)
      % 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)
      % SHOCKBC  Evaluate the residual in the boundary conditions
      res = [ ya(1)+2
         yb(1)  ];
   end
% -----------------------------------------------------------------------

   function jac = shockJac(x,y)
      % 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)
      % 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