常微分方程（十二）：边值问题

初值问题给定起始状态，要求向前推进；边值问题则在两个不同位置提供部分信息，要求找出一条同时满足两端条件的解。措辞上只是微调，后果却大相径庭：边值问题可能有唯一解、无解，甚至无穷多解。它需要一套截然不同的工具箱——本质上是迭代的、全局的，并且与线性代数密不可分。

也正是在这里，常微分方程（ODE）的方法悄然演变为偏微分方程（PDE）的方法。你在本章遇到的离散化、特征值和配点思想，可以直接推广到高维椭圆型 PDE。

从 IVP 到 BVP：小改动，大影响#

相比之下，初值问题会指定 $y(a) = \alpha$ 和 $y'(a) = \alpha'$ 。两者提供的数据量相同（均为两个标量），但数据的分布方式却至关重要。

对于初值问题，在温和的 Lipschitz 条件下，Picard 定理保证了解的存在性与唯一性——解流必然存在。但对于边值问题，这一切都无法保证。

教科书级的经典例子#

考虑方程 $y'' + y = 0$ ，搭配不同的边界条件：

同一个方程，却因边界条件不同而命运迥异。原因在于 $\sin x$ 恰好在 $x = \pi$ 处为零——这意味着对应的齐次问题存在非平凡解。根据 Fredholm 二择一定理，这正是非齐次问题要么超定、要么欠定的临界情形。因此，边值问题的解是否存在、是否唯一，不仅取决于方程本身，更依赖于边界数据。

三类边界条件#

• Dirichlet 条件（指定函数值）：$y(a) = \alpha,\; y(b) = \beta$ 。例如两端温度被固定，是最常见的类型。
• Neumann 条件（指定导数值）：$y'(a) = \alpha,\; y'(b) = \beta$ 。对应热流给定，或梁的一端自由。
• Robin 条件（函数与导数的线性组合）：$\alpha_1 y(a) + \alpha_2 y'(a) = \gamma$ 。例如牛顿冷却定律或部分吸收边界。

实际应用中常混合使用（如一端 Dirichlet，另一端 Neumann）。我们后续发展的数值方法只需稍作调整，即可处理这三类边界条件。

打靶法#

我们的目标是找到 $s$ 使得 $F(s) = 0$ 。这本质上是一个一维求根问题，可用二分法、割线法、牛顿法或 brentq 等方法求解。

这一几何图像恰如炮兵校准火炮：发射不同仰角的炮弹，不断调整直至命中目标——“打靶”之名由此而来。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq

def shooting_method(f, a, b, alpha, beta, s_low=-10, s_high=10):
    """打靶法解 y'' = f(x, y, y'), y(a)=alpha, y(b)=beta。"""
    def residual(s):
        sol = solve_ivp(lambda x, Y: [Y[1], f(x, Y[0], Y[1])],
                        [a, b], [alpha, s], dense_output=True,
                        rtol=1e-9, atol=1e-12)
        return sol.sol(b)[0] - beta

    # 必要时扩张区间
    while residual(s_low) * residual(s_high) > 0:
        s_low, s_high = s_low * 2, s_high * 2

    s_opt = brentq(residual, s_low, s_high, xtol=1e-10)
    sol = solve_ivp(lambda x, Y: [Y[1], f(x, Y[0], Y[1])],
                    [a, b], [alpha, s_opt], dense_output=True,
                    rtol=1e-9, atol=1e-12)
    x_out = np.linspace(a, b, 200)
    return x_out, sol.sol(x_out)[0]

打靶法何时有效，何时失效？#

当底层的初值问题良态时，打靶法表现极佳：$s$ 的微小变化只会引起 $y(b)$ 的微小变化，求根目标清晰明确。然而，若初值问题会指数级放大误差（这在奇异摄动附近或含增长模的刚性问题中很常见），打靶法就会彻底失效。例如，尝试在区间 $[0, 10]$ 上用打靶法求解 $y'' - 100y = 0$ ，边界条件为 $y(0) = 0, y(10) = 1$ ：由于 $\sinh$ 模式的主导作用，残差 $F(s)$ 的取值范围高达 $\sim 10^{43}$ ，使得精确求根几乎不可能。

标准的应对策略是多重打靶法：将区间 $[a, b]$ 划分为 $M$ 个子段，在每段上独立进行打靶（此时每段两端均为未知量），并在相邻子段之间添加匹配条件。虽然最终得到的非线性方程组规模更大，但每一段的积分路径变短，从而保持了良好的数值稳定性。

有限差分法#

将这些近似代入 BVP 在每个内点 $i = 1,\ldots,N-1$ 处的方程，就得到一个包含 $N-1$ 个方程、$N-1$ 个未知数 $y_1, \ldots, y_{N-1}$ 的代数系统（边界值 $y_0 = \alpha$ 和 $y_N = \beta$ 已被纳入右端项）。

对于线性 BVP $y'' + p(x) y' + q(x) y = r(x)$ ，所得系统矩阵是三对角的，仅主对角线及其上下两条对角线非零。求解一个 $n \times n$ 的三对角系统仅需 $\mathcal{O}(n)$ 时间，远优于通用高斯消元法的 $\mathcal{O}(n^3)$ 。

对于非线性 BVP，离散后得到的是非线性代数方程组，通常用牛顿迭代法求解。其雅可比矩阵继承了相同的稀疏结构，因此每次牛顿迭代的成本仍为 $\mathcal{O}(N)$ 。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import numpy as np
from scipy.linalg import solve_banded

def finite_difference_linear(p, q, r, a, b, alpha, beta, N):
    """中心差分求解 y'' + p(x)y' + q(x)y = r(x), y(a)=alpha, y(b)=beta。"""
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1
    pi = np.array([p(xi) for xi in x[1:-1]])
    qi = np.array([q(xi) for xi in x[1:-1]])
    ri = np.array([r(xi) for xi in x[1:-1]])

    main = -2 / h**2 + qi
    upper = 1 / h**2 + pi[:-1] / (2 * h)
    lower = 1 / h**2 - pi[1:]  / (2 * h)

    rhs = ri.copy()
    rhs[0]  -= (1 / h**2 - pi[0]  / (2 * h)) * alpha
    rhs[-1] -= (1 / h**2 + pi[-1] / (2 * h)) * beta

    ab = np.zeros((3, n))
    ab[0, 1:] = upper
    ab[1, :]  = main
    ab[2, :-1] = lower

    y_int = solve_banded((1, 1), ab, rhs)
    return x, np.concatenate([[alpha], y_int, [beta]])

收敛性#

中心差分法是二阶精度的：截断误差为 $\mathcal{O}(h^2)$ ，因此步长 $h$ 减半时，误差大约降至四分之一。若要获得四阶精度，可采用紧致（Pade）有限差分或 Numerov 方法（专为 $y'' = f(x, y)$ 设计）；若追求谱精度，则可使用 Chebyshev 配点法——但代价是矩阵不再稀疏。

特征值问题#

直接求解可得特征值 $\lambda_n = n^2$ ，对应的特征函数为 $\sin(nx)$ ，其中 $n = 1, 2, 3, \ldots$ 。

经有限差分离散后，微分算子 $-y''$ 被替换为一个对称三对角矩阵。该矩阵的特征值近似于理论值 $\{\lambda_n\}$ ：对低阶模态（长波长）精度很高，但随着模态频率升高（振荡尺度接近或小于 $2h$ ），误差逐渐增大，因为网格无法分辨更精细的振荡。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

from scipy.linalg import eigh_tridiagonal
import numpy as np

N = 400
h = np.pi / N
n = N - 1
main = (2 / h**2) * np.ones(n)
off  = (-1 / h**2) * np.ones(n - 1)
eigvals, eigvecs = eigh_tridiagonal(main, off, select='i',
                                    select_range=(0, 9))
print(eigvals.round(4))    # [1.0000  4.0000  9.0000 16.0000 25.0000 ...]

Sturm-Liouville 理论#

其中 $p(x) > 0$ 、$w(x) > 0$ 。这就是著名的 Sturm-Liouville 形式，也是连接常微分方程与数学物理的关键桥梁。

正则 Sturm-Liouville 问题（定义在有限区间上）具有以下性质：

1. 实数、单重、有序的特征值：$\lambda_1 < \lambda_2 < \lambda_3 < \cdots \to \infty$ 。
2. 正交性：当 $m \ne n$ 时，$\int_a^b w(x) y_m(x) y_n(x)\,dx = 0$ 。
3. 完备性：任何满足边界条件的足够光滑函数 $f(x)$ 都可展开为 $f(x) = \sum_n c_n y_n(x)$ ，其中系数 $c_n = \frac{\int_a^b w(x) f(x) y_n(x)\,dx}{\int_a^b w(x) y_n^2(x)\,dx}$ 。
4. 振荡定理：第 $n$ 个特征函数 $y_n(x)$ 在开区间 $(a, b)$ 内恰好有 $n-1$ 个零点。

这些性质正是 Fourier 级数、贝塞尔函数展开、勒让德级数以及量子力学束缚态理论得以成立的基础。它们构成了你未来所有“特征函数展开”工作的正确思维框架。

实例：量子谐振子#

这是一个 Sturm-Liouville 问题，其中 $p \equiv 1$ 、$q(x) = x^2$ 、$w \equiv 1$ ，定义在无界域上。其精确特征值为 $E_n = 2n + 1$ （$n = 0, 1, 2, \ldots$ ），特征函数为厄米特函数。

数值求解时，我们将定义域截断为 $[-L, L]$ （取足够大的 $L$ ，使得波函数在边界处已近似为零），采用中心差分离散，并计算所得矩阵的最低几个特征值。

这一“离散 → 求特征值 → 按能量绘图”的流程，基本上就是处理一维量子束缚态问题的完整工具链。同样的方法也适用于求解振动弦、鼓面或分子的简正模。

配点法：scipy.integrate.solve_bvp#

现代生产级的 BVP 求解器通常采用配点法，而非打靶法或有限差分法。其核心思想是：假设解在网格上是分段多项式，并要求它在每个子区间内的若干“配点”处精确满足微分方程。SciPy 中的 solve_bvp 实现了四阶 Lobatto IIIa 格式，并支持自适应网格加密。

用户需提供：

• fun(x, y)：以向量化形式给出的一阶方程组右端函数；
• bc(ya, yb)：边界条件的残差函数，返回长度等于系统维数的数组；
• x 和 y_init：初始网格和初始猜测解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

from scipy.integrate import solve_bvp
import numpy as np

def fun(x, y):
    return np.vstack([y[1], -2 * np.exp(y[0])])

def bc(ya, yb):
    return np.array([ya[0], yb[0]])      # 两端为零

x_init = np.linspace(0, 1, 10)
y_init = np.zeros((2, x_init.size))
y_init[0] = 0.5 * np.sin(np.pi * x_init)        # 凸起形猜测

sol = solve_bvp(fun, bc, x_init, y_init, tol=1e-10)
assert sol.success

solve_bvp 的优势包括：

• 自适应网格加密：自动在解变化剧烈的区域增加节点；
• 原生支持非线性 BVP：通过牛顿迭代直接求解配点残差系统；
• 连续可微的输出：结果 sol.sol 是一个样条函数，可在任意点求值。

其弱点在于：作为全局方法，它依赖于一个合理的初始猜测；若初始猜测过于离谱，牛顿迭代可能收敛到错误的解分支，甚至完全不收敛。

一题多解：同一问题的三种方法#

在低参数分支上存在两个解（一个小解和一个大解），我们聚焦于小解。

这种多重验证令人安心。当两种独立的 BVP 方法在超出容差范围的情况下仍不一致时，必定至少有一种方法出错了——通常是打靶法遭遇了敏感性壁垒，或是 solve_bvp 因糟糕的初始猜测而陷入错误的解分支。

通向何方：应用#

• 梁的挠度分析：欧拉-伯努利梁方程 $EI\,y^{(4)} = w(x)$ ，结合固定、自由或简支等边界条件。其特征值对应梁的固有振动频率，特征函数即为结构动力学中的模态形状。
• 稳态热传导：$(k(x) T'(x))' + Q(x) = 0$ ，搭配 Dirichlet（指定温度）或 Neumann（指定热流）边界条件。这是典型的 Sturm-Liouville 问题；通过分离变量法，其特征函数展开可求解含时热传导方程。
• 量子束缚态：如前所述。无论是有限深势阱、周期势还是类氢原子的库仑势，数值求解流程完全相同，只需更换势能函数 $V(x)$ 。
• 经典非线性 BVP：Thomas-Fermi、Bratu、Falkner-Skan、Blasius、Painlevé 等问题均有深厚历史背景。只要提供合理的初始猜测，solve_bvp 基本能全部搞定。
• 两点最优控制：庞特里亚金极大值原理将最优控制问题转化为状态-协态系统的 BVP。多重打靶法在此领域是主力工具。

有限差分与配点的思想还可自然推广到二维、三维的椭圆型 PDE（如 Laplace、Poisson、Helmholtz 方程）。此时矩阵虽不再三对角，但仍保持稀疏性；多网格法、Krylov 子空间等迭代求解器可高效处理这类线性代数问题。你在本章学到的一切，都能无缝扩展到更高维度。

方法选择速查表#

实用建议：面对新的 BVP，优先尝试 solve_bvp。若失败，问题通常出在初始猜测上，而非方法本身。

总结#

• 为什么边值问题（BVP）在本质上比初值问题（IVP）更难，并通过一个实例展示存在性与唯一性如何双双失效
• 三种标准边界条件类型：Dirichlet、Neumann 和 Robin
• 打靶法：将 BVP 转化为对缺失初值的根查找问题
• 有限差分法：将 BVP 转化为稀疏线性方程组
• Sturm-Liouville 理论：如何将 BVP 转为矩阵特征值问题（这是通往量子力学与振动分析的桥梁）
• 初探 配点法，通过 scipy.integrate.solve_bvp 这一标准 Python 工具实现

前置知识：需掌握 第 11 章 中介绍的数值方法。

练习题#

1. 存在性与唯一性实验：使用 solve_bvp 尝试求解 BVP $y'' + y = 0$ ，边界条件为 $y(0) = 0$ 、$y(\pi) = 1$ 。记录程序行为，并用 Fredholm 二择一定理解释结果。
2. 有限差分收敛性验证：在 $N \in \{20, 40, 80, 160, 320\}$ 下求解 $y'' = -\pi^2 \sin(\pi x)$ ，$y(0) = y(1) = 0$ 。在双对数坐标下绘制最大误差与步长 $h$ 的关系，并验证斜率为 2。
3. 量子势阱问题：在区间 $[-10, 10]$ 上计算 $-\psi'' + V(x)\psi = E\psi$ 的最低五个特征值与特征函数，其中 $V(x) = -10\,e^{-x^2/2}$ 。该高斯势阱支持多少个束缚态（即 $E < 0$ 的解）？
4. 多重打靶实践：先用单打靶法求解 $y'' = 100\,y$ ，$y(0) = 0$ ，$y(1) = 1$ （应失败或精度极差）。然后将 $[0, 1]$ 分为 5 段，实施多重打靶，并验证结果是否与解析解 $y = \sinh(10x)/\sinh(10)$ 一致。
5. 非线性多解探索：使用 solve_bvp 求解 Bratu 问题 $y'' + \lambda e^y = 0$ ，$y(0) = y(1) = 0$ 在 $\lambda = 2$ 时的两个解分支。（提示：尝试两种不同幅度的“凸起”作为初始猜测。）

下一步#

下一章迈出 ODE 的世界，进入偏微分方程（PDE）。当问题依赖多个独立变量（空间+时间），ODE 不够用了。下一章作为衔接：从波动方程、热方程、Laplace 方程三大经典 PDE 开始，介绍它们的物理意义、典型解法以及与 ODE 的本质区别。

参考文献#

• Ascher, U. M., Mattheij, R. M. M., & Russell, R. D. (1995). Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM Classics.
• Keller, H. B. (1976). Numerical Solution of Two Point Boundary Value Problems. SIAM.
• Boyd, J. P. (2001). Chebyshev and Fourier Spectral Methods (2nd ed.). Dover.
• LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM.
• SciPy documentation: scipy.integrate.solve_bvp.

上一章: Chapter 11: Numerical Methods

下一章: Chapter 13: Introduction to Partial Differential Equations

本文是常微分方程系列（共 18 篇）的第 12 篇。