bvpset matlab,求解边界值问题 - 四阶方法

编写方程代码

使用代换法 y1=y 和 y2=y′，您可以将 ODE 重写为一阶方程组

y1′=y2，

y2′=-2xy2-1x4y1。

对应的函数是

function dydx = bvpfcn(x,y)

dydx = [y(2)

-2*y(2)/x - y(1)/x^4];

end

注意：所有函数都作为局部函数包含在示例末尾。

编写边界条件代码

边界条件函数要求边界条件的形式为 g(y(a),y(b))=0。在此形式中，边界条件是

y(13π)=0，

y(1)-sin(1)=0。

对应的函数是

function res = bcfcn(ya,yb)

res = [ya(1)

yb(1)-sin(1)];

end

设置选项

使用 bvpset 打开求解器统计信息的显示，并指定宽松误差容限以突出求解器之间误差控制的差异。此外，为了提高效率，请指定解析 Jacobian 矩阵

J=∂fi∂y=[∂f1∂y1∂f1∂y2∂f2∂y1∂f2∂y2]=[01-1x4-2x]。

返回 Jacobian 矩阵值的对应函数为

function dfdy = jac(x,y)

dfdy = [0 1

-1/x^4 -2/x];

end

opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');

创建初始估计值

使用 bvpinit 创建解的初始估计值。指定一个常量函数作为初始估计值，初始网格包含区间 [1/3π,1] 中的 10 个点。

xmesh = linspace(1/(3*pi), 1, 10);

solinit = bvpinit(xmesh, [1; 1]);

求解方程

用 bvp4c 和 bvp5c 分别求解方程。

sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 9 points.

The maximum residual is 9.794e-02.

There were 157 calls to the ODE function.

There were 28 calls to the BC function.

sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 11 points.

The maximum error is 6.742e-02.

There were 244 calls to the ODE function.

There were 29 calls to the BC function.

绘制结果

绘制 y1 的两次计算结果和解析解以进行比较。解析解是

y1=sin(1x)，

y2=-1x2cos(1x)。

xplot = linspace(1/(3*pi),1,200);

yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2];

plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo')

title('Comparison of BVP Solvers with Crude Error Tolerance')

legend('True','BVP4C','BVP5C')

xlabel('x')

ylabel('solution y')

绘图验证了 bvp5c 直接控制计算中的真误差，而 bvp4c 仅间接控制计算中的真误差。在更严格的误差容限下，求解器之间的这种差异没有这么明显。

局部函数

此处列出了 BVP 求解器用于求解问题的局部函数。

function dydx = bvpfcn(x,y) % equation to solve

dydx = [y(2)

-2*y(2)/x - y(1)/x^4];

end

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

function res = bcfcn(ya,yb) % boundary conditions

res = [ya(1)

yb(1)-sin(1)];

end

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

function dfdy = jac(x,y) % analytical jacobian for f

dfdy = [0 1

-1/x^4 -2/x];

end

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