求解边界值问题

在边界值问题 (BVP) 中，目标是求常微分方程 (ODE) 的解，该解还需满足某些指定的边界条件。边界条件指定积分区间中两个或多个位置处的解的值之间的关系。在最简单的情形中，边界条件适用于区间的开始和结束（即边界）。

MATLAB^® BVP 求解器 bvp4c 和 bvp5c 用于处理以下形式的 ODE 方程组：

$y' = f (x, y)$

其中：

x 是自变量
y 是因变量
y' 表示 y 关于 x 的导数，也写为 dy/dx

边界条件

在两点 BVP 的最简单情形中，ODE 的解在区间 [a, b] 中求得，并且必须满足边界条件

$g (y (a), y (b)) = 0 .$

要指定给定 BVP 的边界条件，您必须：

编写 res = bcfun(ya,yb) 形式的函数，如果涉及未知参数，则使用 res = bcfun(ya,yb,p) 形式。将此函数作为第二个输入参数提供给求解器。该函数返回 res，这是在边界点处的解的残差值。例如，如果 y(a) = 1 且 y(b) = 0，则边界条件函数为
```
function res = bcfun(ya,yb)
res = [ya(1)-1
       yb(1)];
end
```
在解的初始估计值中，网格中的第一个和最后一个点指定强制执行边界条件的点。对于上述边界条件，您可以指定 bvpinit(linspace(a,b,5),yinit) 来对 a 和 b 强制执行边界条件。

MATLAB 中的 BVP 求解器还适用于包含下列各项的其他类型问题：

未知参数 p
解具有奇异性
多点条件（将积分区间分成若干区域的内边界）

在多点边界条件的情形下，边界条件应用于积分区间中两个以上的点。例如，可能要求在区间的开始处、中间处和结束处的解为零。有关如何指定多点边界条件的详细信息，请参阅 bvpinit。

解的初始估计值

与初始值问题不同，边界值问题的解可以是：

无解
有限个解
无限多个解

求解 BVP 过程的重要部分是提供所需解的估计值。估计值的准确与否对求解器性能甚至是能否成功计算来说至关重要。

使用 bvpinit 函数为解的初始估计值创建结构体。求解器 bvp4c 和 bvp5c 接受此结构体作为第三个输入参数。

为解创建良好的初始估计值更像是一门艺术，而不是一门科学。不过，仍有一些常规指导原则，包括：

确保初始估计值满足边界条件，因为解也需要满足这些边界条件。如果问题包含未知参数，则参数的初始估计值也应该满足边界条件。
尝试将关于实际问题或预期解的信息尽可能多地纳入初始估计值。例如，如果解应该振荡或者有一定数量的符号变化，则初始估计值也应该如此。
考虑网格点的位置（解的初始估计值的 x 坐标）。BVP 求解器会在求解过程中调整这些点，因此您不需要指定太多网格点。最佳做法是在解会快速变化的位置附近指定几个网格点。
如果在较小区间内有一个已知的、更简单的解，则将它用作较大区间内的初始估计值。通常，您可以将一个问题作为一系列相对简单的问题来求解，这种做法称为延拓。使用延拓时，可以使用一个问题的解作为求解下一个问题的初始估计值，从而将一系列简单问题连接起来。

查找未知参数

通常，BVP 会涉及需要同时求解的未知 p 参数。ODE 和边界条件变为

$\begin{array}{l} y' = f (x, y, p) \\ g (y (a), y (b), p) = 0 \end{array}$

在此情况下，边界条件必须足以确定参数 p 的值。

您必须为求解器提供任何未知参数的初始估计值。当调用 bvpinit 来创建结构体 solinit 时，请在第三个输入参数 parameters 中指定初始估计值向量。

solinit = bvpinit(x,v,parameters)

此外，用于编写 ODE 方程和边界条件代码的函数 odefun 和 bcfun 都必须具有第三个参数。

dydx = odefun(x,y,parameters)
res = bcfun(ya,yb,parameters)

求解微分方程时，求解器会调整未知参数的值以满足边界条件。求解器会在 sol.parameters 中返回未知参数的最终值。

奇异 BVP

bvp4c 和 bvp5c 可对以下形式的一类奇异 BVP 求解

$\begin{matrix} y^{'} = \frac{1}{x} S y + f (x, y), \\ 0 = g (y (0), y (b)) . \end{matrix}$

该求解器还可以接受以下形式的问题的未知参数

$\begin{matrix} y^{'} = \frac{1}{x} S y + f (x, y, p), \\ 0 = g (y (0), y (b), p) . \end{matrix}$

奇异问题必须位于 [0,b] 区间上，且 b > 0。使用 bvpset 将常量矩阵 S 作为 'SingularTerm' 选项的值传递给求解器。x = 0 时的边界条件必须与平滑解 Sy(0) = 0 的必要条件一致。解的初始估计值也应该满足此条件。

当您求解奇异 BVP 时，求解器要求您的函数 odefun(x,y) 只返回方程中 f(x, y) 项的值。涉及 S 的项由求解器使用 'SingularTerm' 选项单独处理。

BVP 求解器的选择

MATLAB 提供了求解器 bvp4c 和 bvp5c 来求解 BVP。在大多数情况下，您可以互换使用这些求解器。这些求解器之间的主要区别在于 bvp4c 实现四阶公式，而 bvp5c 实现五阶公式。

除两个求解器之间的误差容限含义以外，bvp5c 函数的使用方法与 bvp4c 完全类似。如果 S(x) 近似于解 y(x)，则 bvp4c 控制残差 |S′(x) – f(x,S(x))|。这种方法间接控制真误差 |y(x) – S(x)|。使用 bvp5c 直接控制真误差。

求解器描述

求解器	描述
`bvp4c`	`bvp4c` 是一个有限差分代码，此代码实现 3 阶段 Lobatto IIIa 公式。这是配置公式，并且配置多项式会提供在整个积分区间中处于四阶精度的 C¹ 连续解。网格选择和误差控制均基于连续解的残差。配置方法使用点网格将积分区间分为子区间。通过对源于边界条件以及所有子区间上配置条件的线性代数方程全局组求解，求解器会确定数值解。然后，求解器会估计每个子区间上数值解的误差。如果解不满足容差标准，则求解器会调整网格并重复计算过程。您必须提供初始网格的点以及网格点处解的初始近似估计。
`bvp5c`	`bvp5c` 是一个有限差分代码，此代码实现 4 阶段 Lobatto IIIa 公式。这是配置公式，并且配置多项式会提供在整个 [a,b] 中具有一致五阶精度的 C¹ 连续解。该公式作为隐式 Runge-Kutta 公式实现。`bvp5c` 直接对代数方程求解；然而，`bvp4c` 使用解析压缩法。`bvp4c` 直接处理未知参数；而 `bvp5c` 使用未知参数的平凡微分方程扩充方程组。

bvp4c

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

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

bvp5c

bvp5c 是一个有限差分代码，此代码实现 4 阶段 Lobatto IIIa 公式。这是配置公式，并且配置多项式会提供在整个 [a,b] 中具有一致五阶精度的 C¹ 连续解。该公式作为隐式 Runge-Kutta 公式实现。bvp5c 直接对代数方程求解；然而，bvp4c 使用解析压缩法。bvp4c 直接处理未知参数；而 bvp5c 使用未知参数的平凡微分方程扩充方程组。

解的计算

在 bvp4c 和 bvp5c 中实现的配置方法在积分区间 [a,b] 上产生 C¹ 连续解。可以使用求解器返回的辅助函数 deval 和结构体 sol，在 [a,b] 中的任意点计算 S(x) 近似解。例如，要在网格点 xint 上计算解 sol，请使用以下命令

Sxint = deval(sol,xint)

deval 函数已向量化。对于向量 xint，Sxint 的第 i 列与解 y(xint(i)) 近似。

BVP 示例和文件

我们可以通过几个示例文件来很好地了解如何求解最常见的 BVP 问题。要方便地查看和运行示例，可以使用 Differential Equations Examples App。要运行此 App，请键入

odeexamples

要打开单独的示例文件进行编辑，请键入

edit exampleFileName.m

要运行示例，请键入

exampleFileName

下表包含可用的 BVP 示例文件及其使用的求解器和选项的列表。

示例文件	使用的求解器	指定的选项	描述	示例链接
`emdenbvp`	`bvp4c` 或 `bvp5c`	`'SingularTerm'`	埃姆登方程，奇异 BVP	求解具有奇异项的 BVP
`fsbvp`	`bvp4c` 或 `bvp5c`	—	无限区间上的 Falkner-Skan BVP	使用延拓验证 BVP 一致性
`mat4bvp`	`bvp4c` 或 `bvp5c`	—	马蒂厄方程的第四个特征函数	求解具有未知参数的 BVP
`rcbvp`	`bvp4c` 和 `bvp5c`	`'FJacobian'` `'AbsTol'` `'RelTol'` `'Stats'`	比较 `bvp4c` 和 `bvp5c` 控制的误差的示例	比较 bvp4c 和 bvp5c 求解器 (`bvp4c`) 比较 bvp4c 和 bvp5c 求解器 (`bvp5c`)
`shockbvp`	`bvp4c` 或 `bvp5c`	`'FJacobian'` `'BCJacobian'` `'Vectorized'`	冲击层位于 x = 0 的解	使用延拓求解 BVP 问题
`twobvp`	`bvp4c`	—	只有两个解的 BVP	对具有两个解的 BVP 求解
`threebvp`	`bvp4c` 或 `bvp5c`	—	三点边界值问题	求解具有多边界条件的 BVP

参考

[1] Ascher, U., R. Mattheij, and R. Russell. “Numerical Solution of Boundary Value Problems for Ordinary Differential Equations.” Philadelphia, PA: SIAM, 1995, p. 372.

[2] 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.

[3] 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.

[4] Shampine, L.F., and J. Kierzenka. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.