bvp4c

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

语法

sol = bvp4c(odefun,bcfun,solinit)

sol = bvp4c(odefun,bcfun,solinit,options)

说明

sol = bvp4c(odefun,bcfun,solinit) 使用 bcfun 描述的边界条件和初始解估计值 solinit 对 odefun 指定的形式为 y′ = f(x,y) 的微分方程组进行积分。使用 bvpinit 函数创建初始估计值 solinit，该函数还定义了要强制应用 bcfun 中边界条件的点。

示例

sol = bvp4c(odefun,bcfun,solinit,options) 还使用由 options（使用 bvpset 函数创建的参量）定义的积分设置。例如，使用 AbsTol 和 RelTol 选项指定绝对误差容限和相对误差容限，或者使用 FJacobian 选项提供 odefun 的解析偏导数。

示例

全部折叠

求解二阶 BVP

打开实时脚本

使用函数在 MATLAB® 中求解二阶 BVP。对于此示例，请使用二阶方程

${y^{'}}^{'} + y = 0$ .

该方程在区间 $[0, π / 2]$ 上定义并受限于边界条件

$y (0) = 0$ ,

$y (π / 2) = 2$ .

要在 MATLAB 中求解此方程，您需要编写一个将方程表示为一阶方程组的函数、一个边界条件函数和一个初始估计值函数。然后 BVP 求解器使用这三个输入来求解方程。

编写方程代码

编写一个用于编写方程代码的函数。使用代换法 $y_{1} = y$ 和 $y_{2} = y^{'}$ 将方程重写为一阶方程组。

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

${y_{2}}^{'} = - y_{1}$ .

对应的函数是

function dydx = bvpfcn(x,y)
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end

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

编写边界条件代码

编写一个函数，以 $g (y (a), y (b)) = 0$ 形式编写边界条件代码。在此形式中，边界条件是

$y (0) = 0$ ,

$y (π / 2) - 2 = 0$ .

对应的函数是

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-2];
end

创建初始估计值

使用 bvpinit 函数为方程的解创建初始估计值。由于该方程将 ${y^{'}}^{'}$ 与 $y$ 联系起来，合理的估计是该解涉及三角函数。使用由积分区间的五个点构成的网格。网格中的第一个和最后一个值是求解器应用边界条件的位置。

初始估计值的函数接受 $x$ 作为输入，并返回对 $y_{1}$ 和 $y_{2}$ 值的估计值。函数是

function g = guess(x)
g = [sin(x)
     cos(x)];
end

xmesh = linspace(0,pi/2,5);
solinit = bvpinit(xmesh, @guess);

求解方程

结合使用 bvp4c 与导函数、边界条件函数和初始估计值来求解问题。

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

对解进行绘图

plot(sol.x, sol.y, '-o')

Figure contains an axes object. The axes object contains 2 objects of type line.

局部函数

此处列出了 bvp4c 用于求解方程的局部函数。

function dydx = bvpfcn(x,y) % equation to solve
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end
%--------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-2];
end
%--------------------------------
function g = guess(x) % initial guess for y and y'
g = [sin(x)
     cos(x)];
end
%--------------------------------

比较 `bvp4c` 和 `bvp5c` 求解器

打开实时脚本

用两个不同求解器以宽松误差容限求解 BVP，并比较结果。

假设使用二阶 ODE

${y^{'}}^{'} + \frac{2}{x} y^{'} + \frac{1}{x^{4}} y = 0$ .

该方程在区间 $[\frac{1}{3 π}, 1]$ 上定义并受限于边界条件

$y (\frac{1}{3 π}) = 0$ ,

$y (1) = \sin (1)$ .

要在 MATLAB® 中求解此方程，您需要编写一个将方程表示为一阶方程组的函数，编写一个边界条件函数，设置一些选项值，并创建初始估计值。然后 BVP 求解器使用这四个输入来求解方程。

编写方程代码

使用代换法 $y_{1} = y$ 和 $y_{2} = y^{'}$ ，您可以将 ODE 重写为一阶方程组

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

${y_{2}}^{'} = - \frac{2}{x} y_{2} - \frac{1}{x^{4}} y_{1}$ .

对应的函数是

function dydx = bvpfcn(x,y)
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end

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

编写边界条件代码

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

$y (\frac{1}{3 π}) = 0$ ,

$y (1) - \sin (1) = 0$ .

对应的函数是

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-sin(1)];
end

设置选项

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

$J = \frac{\partial f_{i}}{\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 \\ - \frac{1}{x^{4}} & - \frac{2}{x} \end{matrix}]$ .

返回雅可比矩阵值的对应函数为

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.

绘制结果

绘制 $y_{1}$ 的两次计算结果和解析解以进行比较。解析解是

$y_{1} = \sin (\frac{1}{x})$ ,

$y_{2} = - \frac{1}{x^{2}} \cos (\frac{1}{x})$ .

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')

Figure contains an axes object. The axes object with title Comparison of BVP Solvers with Crude Error Tolerance, xlabel x, ylabel solution y contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent True, BVP4C, BVP5C.

绘图验证了 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
%---------------------------------

输入参数

全部折叠

`odefun` — 要求解的函数
函数句柄

要求解的函数，指定为定义要积分的函数的函数句柄。odefun 和 bcfun 必须接受相同数量的输入参量。

要编写 odefun 代码，请对标量 x 和列向量 y 使用函数签名 dydx = odefun(x,y)。返回值 dydx 是数据类型为 single 或 double 的列向量，该列向量对应于 f(x,y)。odefun 必须同时接受输入参量 x 和 y，即使其中一个参量未在函数中使用也是如此。

例如，要解算 $y' = 5 y - 3$ ，请使用此函数：

function dydx = odefun(x,y)
dydx = 5*y-3;
end

对于方程组，odefun 的输出为向量。向量中的每个元素是一个方程的解。例如，要求解

$\begin{array}{l} y'_{1} = y_{1} + 2 y_{2} \\ y'_{2} = 3 y_{1} + 2 y_{2} \end{array}$

使用函数：

function dydx = odefun(x,y)
dydx = zeros(2,1);
dydx(1) = y(1)+2*y(2);
dydx(2) = 3*y(1)+2*y(2);
end

bvp4c 也可以求解解具有奇异性或具有多点边界条件的问题。

示例: sol = bvp4c(@odefun, @bcfun, solinit)

未知参数

如果正在求解的 BVP 包含未知参数，您可以使用函数签名 dydx = odefun(x,y,p)，其中 p 是参数值向量。当您使用此函数签名时，BVP 求解器从 solinit 中提供的参数值的初始估计值开始计算未知参数的值。

数据类型: function_handle

`bcfun` — 边界条件
函数句柄

边界条件，指定为计算边界条件中残差的函数句柄。odefun 和 bcfun 必须接受相同数量的输入参量。

要编写 bcfun 代码，请对列向量 ya 和 yb 使用函数签名 res = bcfun(ya,yb)。返回值 res 是数据类型为 single 或 double 的列向量，它对应于边界点处边界条件的残差值。

例如，如果 y(a) = 1 且 y(b) = 0，则边界条件函数为

function res = bcfun(ya,yb)
res = [ya(1)-1
       yb(1)];
end

由于 y(a) = 1，ya(1)-1 在 x = a 点的残差值应为 0。同样，由于 y(b) = 0，yb(1) 在 x = b 点的残差值应为 0。

强制应用边界条件的边界点 x = a 和 x = b 在初始估计值结构体 solinit 中定义。对于两点边界值问题，a = solinit.x(1) 且 b = solinit.x(end)。

示例: sol = bvp4c(@odefun, @bcfun, solinit)

未知参数

如果正在求解的 BVP 包含未知参数，您可以使用函数签名 res = bcfun(ya,yb,p)，其中 p 是参数值向量。当您使用此函数签名时，BVP 求解器从 solinit 中提供的参数值的初始估计值开始计算未知参数的值。

数据类型: function_handle

`solinit` — 解的初始估计值
结构体

解的初始估计值，指定为结构体。使用 bvpinit 创建 solinit。

与初始值问题不同，边界值问题可以无解、有有限数量的解或者有无限数量的解。求解 BVP 过程的重要部分是提供所需解的估计值。估计值的准确与否对求解器性能甚至是能否成功计算来说至关重要。有关创建良好初始估计值的一些指导原则，请参阅解的初始估计值。

示例: sol = bvp4c(@odefun, @bcfun, solinit)

数据类型: struct

`options` — options 结构体
结构体

options 结构体。使用 bvpset 函数创建或修改 options 结构体。

示例: options = bvpset('RelTol',1e-5,'Stats','on') 指定相对误差容限为 1e-5，并打开求解器统计信息的显示。

数据类型: struct

输出参量

全部折叠

`sol` — 解结构体
结构体

解结构体。您可以使用点索引来访问 sol 中的字段，例如 sol.field1。解 (sol.x,sol.y) 在初始网格 solinit.x 中定义的积分区间上是连续的，且在其中有一个连续一阶导数。您可以使用 sol 和 deval 函数来计算该区间中其他点的解。

结构体 sol 包含以下字段。

字段	描述
`x`	`bvp4c` 选择的网格。此网格通常包含不同于初始网格 `solinit.x` 的点。
`y`	`sol.x` 网格点处的 y(x) 逼近值。
`yp`	`sol.x` 网格点处的 y′(x) 逼近值。
`parameters`	`solinit.parameters` 中指定的未知参数的最终值。
`solver`	`'bvp4c'`
`stats`	与解相关的计算成本统计信息：网格点的数量、残差以及调用 `odefun` 和 `bcfun` 的次数。

详细信息

全部折叠

多点边界值问题

对于多点边界值问题，会在积分区间内的几个点上强制应用边界条件。

bvp4c 可以求解多点边界值问题（其中 a = a₀ < a₁ < a₂ < ...< a_n = b 是区间 [a,b] 中的边界点）。点 a₁,a₂,...,a_n−1 表示将 [a,b] 分为多个区域的界点。bvp4c 使用从 1 开始的索引按从左到右的顺序（从 a 到 b）枚举区域。在区域 k [a_k−1,a_k] 中，bvp4c 按如下所示计算导数

yp = odefun(x,y,k)

在边界条件函数 bcfun(yleft,yright) 中，yleft(:,k) 是 [a_k−1,a_k] 的左边界的解。同样，yright(:,k) 是区域 k 右边界的解。具体来说是 yleft(:,1) = y(a) 和 yright(:,end) = y(b)。

当您使用 bvpinit 创建初始估计值时，对于每个交界点，请在 xinit 中使用双重项。有关详细信息，请参阅 bvpinit 的参考页。

如果 yinit 为函数，bvpinit 调用 y = yinit(x,k) 以获取区域 k 中的 x 处的初始估计解。在 bpv4c 返回的解结构体 sol 中，sol.x 为每个交界点指定了双重项。sol.y 的对应列分别包含界点的左解和右解。

有关求解三点边界值问题的示例，请参阅求解具有多边界条件的 BVP。

奇异边界值问题

bvp4c 对一类奇异边界值问题求解，其中包括具有未知参数 p 的形式如下的问题

$\begin{array}{l} y' = S \frac{y}{x} + f (x, y, p), \\ 0 = b c (y (0), y (b), p) . \end{array}$

需要区间为 [0, b] 并且 b > 0。当计算偏微分方程 (PDE) 因柱状或球面对称性而产生的 ODE 的平滑解时，通常会出现此类问题。对于奇异问题，应将（常量）矩阵 S 指定为 bvpset 的 'SingularTerm' 选项的值，并且 odefun 仅计算 f(x,y,p)。边界条件和初始估计值必须与平滑度 S·y(0) = 0 的必要条件一致。

有关求解奇异边界值问题的示例，请参阅求解具有奇异项的 BVP。

算法

bvp4c 是一个有限差分代码，此代码实现 3 阶段 Lobatto IIIa 公式 [1]、[2]。这是配置公式，并且配置多项式会提供在整个积分区间中处于四阶精度的 C¹ 连续解。网格选择和误差控制均基于连续解的残差。

配置方法使用点网格将积分区间分为子区间。通过对源于边界条件以及所有子区间上配置条件的线性代数方程全局组求解，求解器会确定数值解。然后，求解器会估计每个子区间上数值解的误差。如果解不满足容差标准，则求解器会调整网格并重复计算过程。您必须提供初始网格以及网格点处解的初始近似估计。

参考

[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.

[2] Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.

扩展功能

全部展开

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

此函数完全支持基于线程的环境。有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

版本历史记录

在 R2006a 之前推出

另请参阅

主题

求解边界值问题

bvp4c

语法

说明

示例

求解二阶 BVP

比较 bvp4c 和 bvp5c 求解器

输入参数

odefun — 要求解的函数 函数句柄

未知参数

bcfun — 边界条件 函数句柄

未知参数

solinit — 解的初始估计值 结构体

options — options 结构体 结构体

输出参量

sol — 解结构体 结构体

详细信息

多点边界值问题

奇异边界值问题

算法

参考

扩展功能

基于线程的环境 使用 MATLAB® backgroundPool 在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool 加快代码运行速度。

版本历史记录

另请参阅

主题

比较 `bvp4c` 和 `bvp5c` 求解器

`odefun` — 要求解的函数
函数句柄

`bcfun` — 边界条件
函数句柄

`solinit` — 解的初始估计值
结构体

`options` — options 结构体
结构体

`sol` — 解结构体
结构体

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。