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

初值问题给你一个起始状态，让你向前推进；边值问题在两个不同的点上各给你一部分信息，要求你找出两端都吻合的解。措辞改动微小，后果巨大：边值问题可能有唯一解，也可能完全无解，或者有无穷多解。它们要求一套截然不同的工具——迭代的、全局的、与线性代数深度交织的。

这也是 ODE 方法悄悄变成 PDE 方法的地方。这里学到的离散化、特征值、配点思想，可以直接扩展到高维椭圆型 PDE。

本章要点

为什么 BVP 比 IVP 质上更难——存在性和唯一性都可能失败
三种标准边界条件：Dirichlet、Neumann、Robin
打靶法：把 BVP 化简为缺失初值上的求根问题
有限差分：把 BVP 化简为稀疏线性方程组
Sturm-Liouville 理论：BVP 化为矩阵特征值问题（通往量子力学和振动分析的桥梁）
配点法：通过 scipy.integrate.solve_bvp 这一 Python 标准工具

前置知识：第 11 章的数值方法。

1. 从 IVP 到 BVP：小变化，大后果

经典二阶 BVP：$y'' = f(x, y, y'), \quad y(a) = \alpha, \quad y(b) = \beta.$与 IVP（指定$y(a)=\alpha, y'(a)=\alpha'$）相比，数据量相同（两个标量），但分布至关重要。

IVP 在轻度 Lipschitz 条件下由 Picard 定理保证存在唯一性——流是存在的，没什么好说的。BVP 没有任何这种保证。

教科书必讲的例子

考虑$y'' + y = 0$及不同边界条件：

边界条件	解
$y(0)=0,\;y(\pi)=0$	无穷多个：$y = A\sin x$，$A$任意
$y(0)=0,\;y(\pi)=1$	无解（$\sin x$在$x=\pi$处恒为零）
$y(0)=0,\;y(1)=1$	唯一：$y = \sin(x)/\sin(1)$

同一方程，三种迥异命运。原因是$\sin x$恰好在$x=\pi$处为零——这意味着齐次问题有非平凡解，根据 Fredholm 二择一原理，非齐次问题此时要么超定要么欠定。对 BVP，存在唯一性取决于边界数据，不仅仅是方程。

三种边界条件

Dirichlet（函数值）：$y(a)=\alpha,\;y(b)=\beta$。两端温度被固定。最常见。
Neumann（导数值）：$y'(a)=\alpha,\;y'(b)=\beta$。规定热流，或梁的自由端。
Robin（线性组合）：$\alpha_1 y(a) + \alpha_2 y'(a) = \gamma$。Newton 冷却定律、部分吸收。

实际中常出现混合条件（一端 Dirichlet 一端 Neumann）。下面发展的数值方法稍作修改即可处理三种条件。

2. 打靶法

最朴素的想法：把 BVP 转成参数化的 IVP。二阶 BVP 已知$y(a)=\alpha$但不知$y'(a)$。猜一个$s$，正向积分$y'' = f(x, y, y'), \quad y(a) = \alpha, \quad y'(a) = s$至$x=b$，看残差$F(s) := y(b; s) - \beta$。要求$F(s)=0$——一维求根问题，二分、割线、Newton、brentq 都行。

左：BVP$y''+y=0,\;y(0)=0,\;y(\pi/2)=1$。五个不同初始斜率给出五条轨迹，只有$s=1$击中目标。右：残差$F(s) = y(\pi/2; s) - 1$。打靶把 BVP 化为找该曲线的零点。

几何画面就是炮兵测距：发射不同仰角的炮弹，调整直到命中目标。这也是该方法的命名由来。

 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]

打靶法什么时候管用，什么时候不管用？

打靶在底层 IVP 良态时极好：$s$的小变化只引起$y(b)$的小变化，求根目标干净。但当 IVP 指数放大误差时打靶就糟糕——奇异摄动附近、含增长模的刚性问题里常见。试试在$[0,10]$上打靶$y'' - 100y = 0$（$y(0)=0, y(10)=1$）：$\sinh$模会污染一切，残差$F(s)$跨越约 $10^{43}$量级，求根精度根本无从谈起。

标准对策是多重打靶：把$[a, b]$切成$M$段，各段独立打靶（两端均为未知），段间加匹配方程。结果是更大的非线性方程组，但每段都短而良态。

3. 有限差分法

不同哲学：先离散算子，再求解。在均匀网格$x_i = a + ih,\;i=0,\ldots,N,\;h = (b-a)/N$上用中心差分$y'(x_i) \approx \frac{y_{i+1} - y_{i-1}}{2h}, \qquad y''(x_i) \approx \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}.$代入 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)$，远低于通用 Gauss 消元的$\mathcal{O}(n^3)$。

左：$y''=-\pi^2\sin(\pi x)$（$y(0)=y(1)=0$）的数值解，精确解为$\sin(\pi x)$。三个高亮点显示局部三点模板。右：$h^2A$结构。只有三条非零对角，其余严格为零。内存与 CPU 开销均为$\mathcal{O}(N)$。

非线性 BVP 离散后得非线性代数系统，用 Newton 迭代求解。Newton-Jacobi 矩阵继承同样的稀疏性，每次 Newton 步仍是$\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 配点——代价是稠密矩阵。

4. 特征值问题

齐次边界条件下含参数$\lambda$的线性 BVP 就是特征值问题。最简单的例子：$-y'' = \lambda y, \quad y(0) = y(\pi) = 0.$直接计算给出特征值$\lambda_n = n^2$与特征函数$\sin(nx)$（$n=1,2,3,\ldots$）。

有限差分离散后算子$-y''$变成对称三对角矩阵。其特征值近似$\{\lambda_n\}$，对低模（长波长）精确，对高模（网格无法分辨小于$2h$的振荡）精度逐渐下降。

左：前五个数值特征模（彩色实线）与精确$\sin(nx)$（黑虚线），错位以便观看。数值特征值精确到小数四位匹配$1, 4, 9, 16, 25$。右：低模严格落在$\lambda_n=n^2$曲线上；系统色散误差随$n$接近网格极限而增大。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
from scipy.linalg import eigh_tridiagonal
import numpy as np

# -y'' = lambda y, y(0) = y(pi) = 0, N 段
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 ...]

5. Sturm-Liouville 理论

一般的二阶自共轭特征值问题：$-\bigl(p(x)\,y'\bigr)' + q(x)\,y = \lambda\,w(x)\,y, \quad y(a) = y(b) = 0,$其中$p>0,w>0$。这就是 Sturm-Liouville 形式——从 ODE 通往数学物理的桥梁。

性质（有限区间上的正则 SL 问题）：

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

正是这些性质使 Fourier 级数、Bessel 展开、Legendre 级数、量子力学束缚态都"能用"。它们是思考一切"特征函数展开"的正确方式。

实例：量子谐振子

量子谐振子的无量纲 Schrödinger 方程：$-\psi''(x) + x^2\,\psi(x) = E\,\psi(x), \quad \psi(\pm\infty) = 0.$这是$p \equiv 1, q(x) = x^2, w \equiv 1$的 SL 问题（无界域）。精确特征值：$E_n = 2n + 1$（$n=0,1,2,\ldots$），特征函数为 Hermite 函数。

数值上把域截断到$[-L, L]$（$L$足够大使$\psi$在边界已几乎为零），中心差分离散，求矩阵的若干最低特征值。

黑曲线为势能$V(x)=x^2$。每个彩色特征函数画在能量$E_n$高度上，呈现经典结构：每升高一模多一个零点；波函数在经典回弹点之间振荡，外侧指数衰减。数值特征值与精确$2n+1$匹配到小数四位。

这一模板——“离散，eigh，按能量画”——基本上就是一维束缚态量子力学的整套工具。同样机制给出振动弦、鼓、分子的简正模。

6. 配点法：scipy.integrate.solve_bvp

现代生产级 BVP 求解器通常用配点法而非打靶或有限差分。思路：假设解是网格上的分段多项式，并要求它在每段内的若干"配点"上严格满足 ODE。SciPy 的 solve_bvp 实现 4 阶 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
16
from scipy.integrate import solve_bvp
import numpy as np

# Bratu 问题：y'' + lambda * exp(y) = 0, y(0) = y(1) = 0, lambda = 2
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 优势：

自适应网格细化——节点放在解变化最快处；
通过 Newton 迭代原生处理非线性 BVP；
结果 sol.sol 是连续可微的样条，可在任意点求值。

弱点：与所有全局方法一样，需要合理的初始猜测；糟糕的猜测会让 Newton 收敛到不同的解分支（或不收敛）。

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

著名的 Bratu 问题$y'' + 2 e^y = 0,\;y(0)=y(1)=0$在低分支上有两个解（小的和较大的），我们关注小的。

左：打靶 + brentq、有限差分 + Newton、solve_bvp 配点——三种不同算法、三条不同数值流水线、一个解。右：两两差异在$10^{-7}$量级，由各方法自身容差决定。

这种冗余令人安心：当两种独立 BVP 方法在容差以外水平不一致时，有一个肯定是错的——通常是打靶撞上敏感性墙，或 solve_bvp 因坏猜测被困在错误分支。

8. 通向何方：应用

梁的挠度。Euler-Bernoulli 梁$EI\,y^{(4)} = w(x)$，配各种固定/自由/简支端组合。特征值给出梁的固有振动频率，特征函数为模态形状（结构动力学的基础）。
稳态热传导。$(k(x)T'(x))' + Q(x) = 0$，配 Dirichlet（指定温度）或 Neumann（指定热流）边界。彻头彻尾的 Sturm-Liouville；特征函数展开通过分离变量得到时间依赖的热方程。
量子束缚态。如上。换有限井、周期势、类氢库仑尾——数值流程不变，只换$V(x)$。
Thomas-Fermi、Bratu、Falkner-Skan、Blasius、Painleve。经典非线性 BVP，历史悠久。给个合理初猜，全部交给 solve_bvp。
两点最优控制。Pontryagin 极大值原理把最优控制化为状态-协态系统的 BVP。多重打靶是这里的工作马。

同样的有限差分与配点思想直接扩展到二维、三维椭圆 PDE（Laplace、Poisson、Helmholtz）。矩阵不再是三对角但仍然稀疏；迭代求解器（multigrid、Krylov）处理线性代数。关于 BVP 学到的一切都可以放大。

9. 方法选择速查

方法	适用	注意
打靶	良态 IVP，简单问题	敏感性 / 刚性爆炸
多重打靶	敏感问题、两点最优控制	设置复杂、非线性系统更大
有限差分	线性问题、特征值问题	阶数受限于二（提阶要额外工作）
谱 / Chebyshev	光滑解、高精度	需要谨慎选基
配点（`solve_bvp`）	通用非线性 BVP	需要合理初猜

实用默认：新 BVP 直接用 solve_bvp，失败往往是初猜的锅，不是方法的锅。

习题

存在唯一性实验。用 solve_bvp 求$y''+y=0,\;y(0)=0,\;y(\pi)=1$。记录会发生什么并用 Fredholm 二择一原理解释。
有限差分收敛性。在$N\in\{20,40,80,160,320\}$下解$y''=-\pi^2\sin(\pi x),\;y(0)=y(1)=0$。log-log 画最大误差对$h$的曲线，验证斜率为 2。
量子势井。在$[-10, 10]$上对$-\psi''+V(x)\psi=E\psi$（$V(x)=-10e^{-x^2/2}$）求最低五个特征值与特征函数。这个高斯井支持几个束缚态（$E<0$）？
多重打靶。用单打靶解$y''=100y,\;y(0)=0,\;y(1)=1$（应失败或精度极差）；再把$[0,1]$分成 5 段做多重打靶，验证结果与$y=\sinh(10x)/\sinh(10)$匹配。
非线性示例。用 solve_bvp 找 Bratu 问题$y''+\lambda e^y=0,\;y(0)=y(1)=0$（$\lambda=2$）的两个分支。提示：试两个不同凸起幅度的初猜。

参考文献

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（第 2 版）. Dover.
LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM.
SciPy 文档：scipy.integrate.solve_bvp。

系列导航


上一章	第 11 章：数值方法
当前	第十二章：边值问题
下一章	第 13 章：偏微分方程引论