MATLAB 数学应用 微分方程 边界值问题 求解边界值问题

在边界值问题 (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)-1yb(1)];end

在解的初始估计值中，网格中的第一个和最后一个点指定强制执行边界条件的点。对于上述边界条件，您可以指定 bvpinit(linspace(a,b,5),yinit) 来对 a 和 b 强制执行边界条件。

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

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

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

解的初始估计值

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

无解有限个解无限多个解

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

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

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

确保初始估计值满足边界条件，因为解也需要满足这些边界条件。如果问题包含未知参数，则参数的初始估计值也应该满足边界条件。

尝试将关于实际问题或预期解的信息尽可能多地纳入初始估计值。例如，如果解应该振荡或者有一定数量的符号变化，则初始估计值也应该如此。

考虑网格点的位置（解的初始估计值的 x 坐标）。BVP 求解器会在求解过程中调整这些点，因此您不需要指定太多网格点。最佳做法是在解会快速变化的位置附近指定几个网格点。

如果在较小区间内有一个已知的、更简单的解，则将它用作较大区间内的初始估计值。通常，您可以将一个问题作为一系列相对简单的问题来求解，这种做法称为延拓。使用延拓时，可以使用一个问题的解作为求解下一个问题的初始估计值，从而将一系列简单问题连接起来。

查找未知参数

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

y′=f(x,y,p)

g(y(a),y(b),p)=0

在此情况下，边界条件必须足以确定参数 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 求解

y′=1xSy+f(x,y)y' = \dfrac{1}{x}Sy + f(x, y)y′=x1​Sy+f(x,y)

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

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

y′=1xSy+f(x,y,p)y' = \dfrac{1}{x}Sy + f(x,y,p)y′=x1​Sy+f(x,y,p)

0=g(y(0),y(b),p).0=g(y(0),y(b),p).0=g(y(0),y(b),p).

奇异问题必须位于 [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 和 bvp5c 中实现的配置方法在积分区间 [a,b] 上产生 C1 连续解。可以使用求解器返回的辅助函数 deval 和结构体 sol，在 [a,b] 中的任意点计算 S(x) 近似解。例如，要在网格点 xint 上计算解 sol，请使用以下命令

Sxint = deval(sol,xint)

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