bvp4c-shockbvp

方程：e*y'' + x*y' = -e*pi^2*cos(pi*x) - pi*x*sin(pi*x)   0<e<<1时

边界条件：y(-1) = -2 and y(1) = 0

开启vectorized

微分方程

function dydx=ode(x,y,e)
%开启vectorized,能够减少运行时间
%和没有开启相比,把相应的y(m)变成y(m,:),*和/换成点乘点即可
pix=pi*x
dydx=[y(2,:) 
    -x.*y(2,:)/e+pi^2*cos(pix)-pix.*sin(pix)/e];

边界条件

function res=bcfun(ya,yb,e)
res=[ya(1)+2
    yb(1)];

微分方程的雅克比矩阵
% [df1/dy(1)  df1/dy(2)
%  df2/dy(1)  df2/dy(2)]
function JJ=shockJac(x,y,e) %e为已知的传递参数
JJ=[0 1
    0 -x/e];

边界条件的雅克比矩阵,两个矩阵
% dBC/dya
% dBC/dyb
function [dBCdya,dBCdyb]=shockBCjac(ya,yb,e)
dBCdya=[1 0
    0 0];
dBCdyb=[0 0
    1 0];

主程序

options = bvpset('FJacobian',@shockJac,'BCJacobian',@shockBCJac,'Vectorized','on');
sol= bvpinit([-1 -0.5 0 0.5 1],[1 0]);
e=1e-4
sol = bvp4c(@ode,@bcfun,sol,options,e);
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');

自己编写的源文件

function example
options = bvpset('FJacobian',@shockJac,'BCJacobian',@shockBCJac,'Vectorized','on');
sol= bvpinit([-1 -0.5 0 0.5 1],[1 0]);
e=1e-4
sol = bvp4c(@ode,@bcfun,sol,options,e);
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');
%微分方程
function dydx=ode(x,y,e)
%开启vectorized
%和没有开启相比,把相应的y(m)变成y(m,:),*和/换成点乘点即可
pix=pi*x
dydx=[y(2,:) 
    -x.*y(2,:)/e+pi^2*cos(pix)-pix.*sin(pix)/e];
   
%边界条件
function res=bcfun(ya,yb,e)
res=[ya(1)+2
    yb(1)];

%微分方程的雅克比矩阵
% [df1/dy(1)  df1/dy(2)
%  df2/dy(1)  df2/dy(2)]
function JJ=shockJac(x,y,e)
JJ=[0 1
    0 -x/e];

%边界条件的雅克比矩阵,两个矩阵
% dBC/dya
% dBC/dyb
function [dBCdya,dBCdyb]=shockBCJac(ya,yb,e)
dBCdya=[1 0
    0 0];
dBCdyb=[0 0
    1 0];
matlab提供的源文件：MATLAB6p5\toolbox\matlab\demos\shockbvp.m