使用延拓求解 BVP 问题

编写方程代码

使用代换法 ${\mathit{y}}_1=\mathit{y}$ 和 ${\mathit{y}}_2={\mathit{y}}^{\prime }$，您可以将方程重写为一阶方程组

${{\mathit{y}}_1}^{\prime }={\mathit{y}}_2$,

${{\mathit{y}}_2}^{\prime }=-\frac{\mathit{x}}{\mathit{e}}{\mathit{y}}^{\prime }-{\pi }^2\cos \left(\pi \mathit{x}\right)-\frac{\pi \mathit{x}}{\mathit{e}}\sin \left(\pi \mathit{x}\right)$.

编写一个函数以使用签名 dydx = shockode(x,y) 编写方程代码，其中：

将函数向量化，以使 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 求解器要求边界条件采用 $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$ 形式。在此形式中，边界条件是：

$\mathit{y}\left(-1\right)+2=0$,

$\mathit{y}\left(1\right)=0$.

编写一个函数以使用签名 res = shockbc(ya,yb) 来编写边界条件代码，其中：

对应的函数是

function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end

编写雅可比矩阵代码

在此问题中，ODE 函数和边界条件的解析雅可比矩阵可以很轻松地计算出来。提供雅可比矩阵使得求解器效率更高，因为求解器不再需要通过有限差分来逼近它们。

对于 ODE 函数，雅可比矩阵为

${\mathit{J}}_{\mathrm{ODE}}=\frac{\partial \mathit{f}}{\partial \mathit{y}}=\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_2}\\ \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& -\frac{\mathit{x}}{\mathit{e}}\end{array}\right]$.

对应的函数是

function jac = shockjac(x,y,e)
jac = [0   1
       0  -x/e];
end

同样，对于边界条件，雅可比矩阵为

${\mathit{J}}_{\mathit{y}\left(\mathit{a}\right)}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$ , ${\mathit{J}}_{\mathit{y}\left(\mathit{b}\right)}=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]$.

对应的函数是

function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end

获取初始估计值

使用常量估计值在包含 $\left[-1,1\right]$ 中的五个点的网格上求解。

sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

求解方程

如果您尝试使用 $\mathit{e}={10}^{-4}$ 直接求解方程，则求解器会由于问题在转变点 $\mathit{x}=0$ 附近处的不良条件而难以求解。在这种情况下，为了获得 $\mathit{e}={10}^{-4}$ 的解，此示例使用了延拓，即对 ${10}^{-2}$、${10}^{-3}$ 和 ${10}^{-4}$ 求解一系列问题。在每次迭代中求解器的输出充当下一次迭代中解的估计值（这就是为什么 bvpinit 的初始估计值的变量是 sol，求解器的输出也命名为 sol）。

由于雅可比矩阵的值取决于 $\mathit{e}$ 的值，因此需要设置循环中的选项，为雅可比矩阵指定 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

对解进行绘图

基于网格 $\mathit{x}$ 和解 $\mathit{y}\left(\mathit{x}\right)$ 绘制 bvp4c 的输出。使用延拓时，求解器能够处理在 $\mathit{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
%-------------------------------------------