图像延拓matlab,使用延拓求解 BVP 问题

编写方程代码

使用代换法 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

%-------------------------------------------