MATLAB 数学应用 微分方程 边界值问题 使用延拓求解BVP问题

本文讲述了如何使用延拓求解难以进行数值求解的边界值问题，延拓实际上是将问题分解成一系列更简单的问题。

对于 0<e≪1，考虑如下微分方程

$e{y}^{\prime {}^{\prime }}+x{y}^{\prime }=-e{\pi }^{2}cos(\pi x)-\pi xsin(\pi x)$。

此问题位于区间 [−1,1] 上，并且需要满足边界条件

y(−1)=−2，

y(1)=0。

当 $e=10^{-4}$ 时，方程的解会在 x=0 附近快速转变，因此难以进行数值求解。此示例使用延拓对 e 的几个值进行迭代处理，直到 $e=10^{-4}$ 。每个中间解都用作下一个问题的初始估计值。

要在 MATLAB 中对此方程组求解，您需要先编写方程组、边界条件和初始估计值的代码，然后再调用边界值问题求解器 bvp4c。您可以将所需的函数作为局部函数包含在文件末尾，也可以将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

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

$y_1^{\prime }=y_2$ ，

$y_2^{\prime }=-\frac{x}{e}{y}^{\prime }-{\pi }^{2}cos(\pi x)-\frac{\pi x}{e}sin(\pi x)$。

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

x 是自变量。

y 是因变量。

dydx(1) 给出 $y_1^{\prime }$ 的方程，dydx(2) 给出 $y_2^{\prime }$ 的方程。

将函数向量化，以使 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 矩阵为

$$J_{OED}=\frac{\partial f}{\partial y}=\left[\begin{array}{cc}\frac{\partial f_1}{\partial y_1}& \frac{\partial f_1}{\partial y_2}\\ \frac{\partial f_2}{\partial y_1}& \frac{\partial f_2}{\partial y_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& -\frac{x}{e}\end{array}\right]$$。

对应的函数是

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

同样，对于边界条件，Jacobian 矩阵为

$J_{y(a)}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$、$J_{y(b)}=\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

获取初始估计值

使用常量估计值在包含 [−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
%-------------------------------------------