编写方程代码
使用代换法 y1=y 和 y2=y′，您可以将方程重写为一阶方程组
y1′=y2，
y2′=-xey′-π2cos(πx)-πxesin(πx)。
编写一个函数以使用签名 dydx = shockode(x,y) 编写方程代码，其中：
x 是自变量。
y 是因变量。
dydx(1) 给出 y1′ 的方程，dydx(2) 给出 y2′ 的方程。
将函数向量化，以使 shockode([x1 x2 …],[y1 y2 …]) 返回 [shockode(x1,y1) shockode(x2,y2) …]。这种方法提高了求解器的性能。
对应的函数是
function dydx = shockode(x,y)
pix = pi*x;
dydx = [y(2,:)
-x/e.*y(2,:) – pi^2*cos(pix) – pix/e.*sin(pix)];
end
注意：所有函数都作为局部函数包含在示例的末尾。
编写边界条件代码
BVP 求解器要求边界条件采用 g(y(a),y(b))=0 形式。在此形式中，边界条件是：
y(-1)+2=0，
y(1)=0。
编写一个函数以使用签名 res = shockbc(ya,yb) 来编写边界条件代码，其中：
ya 是在区间 [a,b] 开始处的边界条件的值。
yb 是在区间 [a,b] 结束处的边界条件的值。
对应的函数是
function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
yb(1)];
end
编写 Jacobian 矩阵代码
在此问题中，ODE 函数和边界条件的解析 Jacobian 矩阵可以很轻松地计算出来。提供 Jacobian 矩阵使得求解器效率更高，因为求解器不再需要通过有限差分来逼近它们。
对于 ODE 函数，Jacobian 矩阵为
JODE=∂f∂y=[∂f1∂y1∂f1∂y2∂f2∂y1∂f2∂y2]=[010-xe]。
对应的函数是
function jac = shockjac(x,y,e)
jac = [0 1
0 -x/e];
end
同样，对于边界条件，Jacobian 矩阵为
Jy(a)=[1000]、Jy(b)=[0010]。
对应的函数是
function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end
获取初始估计值
使用常量估计值在包含 [-1,1] 中的五个点的网格上求解。
sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);
求解方程
如果您尝试使用 e=10-4 直接求解方程，则求解器会由于问题在转变点 x=0 附近处的不良条件而难以求解。在这种情况下，为了获得 e=10-4 的解，此示例使用了延拓，即对 10-2、10-3 和 10-4 求解一系列问题。在每次迭代中求解器的输出充当下一次迭代中解的估计值(这就是为什么 bvpinit 的初始估计值的变量是 sol，求解器的输出也命名为 sol)。
由于 Jacobian 矩阵的值取决于 e 的值，因此需要设置循环中的选项，为 Jacobian 矩阵指定 shockjac 和 shockbcjac 函数。此外，还要启用向量化，因为编写的 shockode 用于处理值向量。
e = 0.1;
for i = 2:4
e = e/10;
options = bvpset(‘FJacobian’,@(x,y) shockjac(x,y,e),’BCJacobian’,@shockbcjac,’Vectorized’,’on’);
sol = bvp4c(@(x,y) shockode(x,y,e),@shockbc, sol, options);
end
对解进行绘图
基于网格 x 和解 y(x) 绘制 bvp4c 的输出。使用延拓时，求解器能够处理在 x=0 处的不连续性。
plot(sol.x,sol.y(1,:),’-o’);
axis([-1 1 -2.2 2.2]);
title([‘There Is a Shock at x = 0 When e =’ sprintf(‘%.e’,e) ‘.’]);
xlabel(‘x’);
ylabel(‘solution y’);
局部函数
此处列出了 BVP 求解器 bvp4c 为计算解而调用的局部函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function dydx = shockode(x,y,e) % equation to solve
pix = pi*x;
dydx = [y(2,:)
-x/e.*y(2,:) – pi^2*cos(pix) – pix/e.*sin(pix)];
end
%——————————————-
function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
yb(1)];
end
%——————————————-
function jac = shockjac(x,y,e) % jacobian of shockode
jac = [0 1
0 -x/e];
end
%——————————————-
function [dBCdya,dBCdyb] = shockbcjac(ya,yb) % jacobian of shockbc
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end
%——————————————-