使用延拓求解 BVP 问题

以下示例说明如何使用延拓求解难以进行数值求解的边界值问题，延拓实际上是将问题分解成一系列更简单的问题。

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

$e {y^{'}}^{'} + {xy}^{'} = - e π^{2} \cos (π x) - π xsin (π 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^{'}$ ，您可以将方程重写为一阶方程组

${y_{1}}^{'} = y_{2}$ ,

${y_{2}}^{'} = - \frac{x}{e} y^{'} - π^{2} \cos (π x) - \frac{π x}{e} \sin (π x)$ .

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

x 是自变量。
y 是因变量。
dydx(1) 给出 ${y_{1}}^{'}$ 的方程，dydx(2) 给出 ${y_{2}}^{'}$ 的方程。

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

编写雅可比矩阵代码

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

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

$J_{ODE} = \frac{\partial f}{\partial y} = [\begin{matrix} \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{matrix}] = [\begin{matrix} 0 & 1 \\ 0 & - \frac{x}{e} \end{matrix}]$ .

对应的函数是

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

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

$J_{y (a)} = [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}]$ , $J_{y (b)} = [\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}]$ .

对应的函数是

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）。

由于雅可比矩阵的值取决于 $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

对解进行绘图

基于网格 $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');

Figure contains an axes object. The axes object with title There Is a Shock at x = 0 When e =1e-04., xlabel x, ylabel solution y contains an object of type line.

局部函数

此处列出了 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
%-------------------------------------------

另请参阅

bvp4c | bvpinit | bvpset

使用延拓求解 BVP 问题

编写方程代码

编写边界条件代码

编写雅可比矩阵代码

获取初始估计值

求解方程

对解进行绘图

局部函数

另请参阅

相关主题

WeChat