常微分方程（十三）：偏微分方程引论

当一个量依赖于不止一个自变量，整个 ODE 世界就分裂为一个远更丰富的世界：偏微分方程（PDE）。 金属棒里的温度同时是位置和时间的函数；振动的弦在空间与时间两个维度中演化；静电势驻留在三维空间中。ODE 的所有技术此时变成"工具"而不是"答案"——分离变量法把一个 PDE 拆成一族 ODE，那族 ODE 的本征值就是算子的谱，叠加原理再把一切重新缝合。

本章不是百科全书，而是一次实战导览：我们聚焦三大经典方程——热方程、波动方程、Laplace 方程——因为任何二阶线性 PDE 在二变量情形下，从一个精确意义上看，都"相似"于三者之一。

本章要点

PDE 与 ODE 的本质差别，以及"三大正则型"为何几乎覆盖一切线性二阶问题
分类定理：判别式 $B^2 - AC$ 决定抛物型 / 双曲型 / 椭圆型
分离变量法的标准流程：PDE -> Sturm-Liouville 本征问题 -> Fourier 级数
d’Alembert 公式与光锥：有限速度传播 vs 形式上的"无限速度"
有限差分法（FTCS、Crank-Nicolson、leapfrog）及其稳定性（CFL 条件、von Neumann 分析）
极值原理 与 Gauss-Seidel 迭代求解 Laplace 方程
适定性 的初步认识：哪些初/边值数据让 PDE 问题"恰好可解"

先修知识：第 11 章数值方法，第 12 章边值问题，以及多元微积分（梯度、Laplace 算子、链式法则）。

从 ODE 到 PDE

ODE 中，未知函数 $y(t)$ 依赖单一变量。每个解由有限多个常数（数目等于阶数）刻画，初值便能挑出唯一一条解。

PDE 要求一个多变量函数。要确定它，需要在整条低维曲面（2D 中的曲线、3D 中的曲面）上给定数据；解空间真正是无限维的。同一个方程，给定一种数据时适定，换一种数据可能立刻不适定——这是 ODE 完全没有的现象。

三大经典方程

方程	公式	类型	物理意义
热方程	$u_t = \alpha\,u_{xx}$	抛物型	不可逆扩散
波动方程	$u_{tt} = c^2\,u_{xx}$	双曲型	时间可逆传播
Laplace 方程	$u_{xx} + u_{yy} = 0$	椭圆型	平衡 / 无时间演化

为什么恰好是这三个？因为任意二变量二阶线性 PDE，经过坐标变换，都可以化为这三种正则型之一。判别式决定属于哪一种。

分类定理

$$A\,u_{xx} + 2B\,u_{xy} + C\,u_{yy} + (\text{低阶项}) = 0,$$

定义判别式 $\Delta = B^2 - AC$：

$\Delta$	类型	正则代表	实特征曲线
$> 0$	双曲型	波动方程	两族
$= 0$	抛物型	热方程	一族（退化）
$< 0$	椭圆型	Laplace 方程	无（复）

实特征曲线的数目，决定扰动如何传播、问题接受什么样的数据、什么样的数值格式才能用。

左上：$(AC, B^2)$ 平面被直线 $B^2 = AC$ 切分——上方是双曲型，下方是椭圆型，线上是抛物型。中上：抛物型基本解——原点处的 $\delta$ 源迅速展开为高斯分布，形式上一瞬间在所有点都非零（形式"无限速度"）。右上：双曲型——局部扰动严格被光锥 $|x|\leq ct$ 限制。下排：三大正则方程并排呈现，每张卡片列出最关键的性质。

物理直觉：

双曲型：池塘里投石子——涟漪以有限速度向外传播；信号到达远处需要时间。
抛物型：触摸钢棒一端——理论上另一端"立刻"开始升温（虽然指数级地微弱）。
椭圆型：稳态——时间已被积出，只剩源与汇之间的平衡。

抛物型的"无限速度"在物理上不严格，但实际上误差衰减极快、不会在可观测时间内显现，所以模型仍然非常成功。

热方程

一段话推导

$$\rho c_p \frac{\partial u}{\partial t} = -\frac{\partial q}{\partial x},$$$$u_t = \alpha\,u_{xx}, \qquad \alpha = \frac{k}{\rho c_p}.$$

扩散系数 $\alpha$ 的量纲是 $\text{长度}^2/\text{时间}$——这是 $k, \rho, c_p$ 唯一自然的组合方式。

分离变量法

$$\frac{T'}{\alpha T} = \frac{X''}{X} = -\lambda \quad (\text{常数}).$$

左边只含 $t$，右边只含 $x$，相等只能是同一个常数。于是

$$X'' + \lambda X = 0, \qquad T' + \alpha\lambda\,T = 0.$$

Dirichlet 边界给出本征值 $\lambda_n = (n\pi/L)^2$，本征函数 $X_n(x) = \sin(n\pi x/L)$，时间因子 $T_n(t) = e^{-\alpha\lambda_n t}$。叠加得到完整的 Fourier 正弦级数解：

$$\boxed{\;u(x, t) = \sum_{n=1}^\infty b_n\,\sin\!\frac{n\pi x}{L}\,e^{-\alpha (n\pi/L)^2 t},\quad b_n = \frac{2}{L}\int_0^L f(x)\sin\!\frac{n\pi x}{L}\,dx.\;}$$

两条结构性教训：

热方程是 Fourier 滤波器。高模态（大 $n$）衰减比低模态快得多——模态 $n$ 的半衰期 $\propto 1/n^2$。初值中的尖锐特征几乎瞬间被磨平，平滑特征则慢慢褪去。
边界条件选择基函数：Dirichlet -> 正弦；Neumann -> 余弦；周期 -> 复指数。本征函数正交，系数即内积。这就是 Laplace 算子的谱分解。

左：$u(x, t)$ 在六个时刻的快照——分段常数初值在不到一秒内被磨平，趋向平衡 $u\equiv 0$。中：同一解的时空热图——明亮结构沿时间轴上升时迅速褪色。右：Fourier 模态在 log 坐标下的衰减；模态 $n=10$ 比 $n=1$ 快上千倍。热方程是一个低通滤波器。

分离变量过程的可视化

把"配方"中的每个原料分别画出来更直观。本征函数是空间驻波；时间因子是衰减包络；它们的张量积是基本解；叠加重构初值。

左上：空间本征模 $X_n(x) = \sin(n\pi x/L)$，自动满足两端为零。右上：对应的时间因子 $T_n(t)$，$n$ 越大衰减越快。中排：基本解 $u_n(x, t) = X_n T_n$ 在三个时刻的形态——短波模态在 $t=2$ 时几乎消失。下排：用越来越多的模态重建一个方波初值。三个模态时形状勉强可见；三十个模态时除了跳跃点附近以外几乎处处吻合，跳跃点处的过冲就是著名的 Gibbs 现象。

热方程的有限差分

离散化 $x_j = j\,\Delta x$、$t^n = n\,\Delta t$。三种主力格式，每一种都是一行更新公式加一段稳定性故事：

$$u_j^{n+1} = u_j^n + r\,(u_{j+1}^n - 2u_j^n + u_{j-1}^n), \qquad r = \frac{\alpha\,\Delta t}{\Delta x^2}.$$

von Neumann 分析（代入 $u_j^n = \xi^n e^{ikj\Delta x}$，要求 $|\xi|\leq 1$）给出 $r \leq 1/2$。$\Delta x$ 减半要求 $\Delta t$ 减为四分之一——精度便宜，稳定性昂贵。

BTCS（后向时间，中心空间）：隐式，$O(\Delta t, \Delta x^2)$，无条件稳定——任何 $r > 0$ 都满足 $|\xi|\leq 1$。每步求解一个三对角系统。

Crank-Nicolson：FTCS 与 BTCS 的平均；$O(\Delta t^2, \Delta x^2)$，无条件稳定，是线性热方程的"金标准"。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
from scipy.linalg import solve_banded

def crank_nicolson_heat(f, alpha, L, T_end, Nx=101, Nt=400):
    """Crank-Nicolson 求解 u_t = alpha u_xx，Dirichlet 边界 u(0)=u(L)=0。"""
    x  = np.linspace(0, L, Nx); dx = x[1] - x[0]
    dt = T_end / Nt
    r  = alpha * dt / dx**2
    u  = f(x); u[0] = 0; u[-1] = 0
    n  = Nx - 2
    ab = np.zeros((3, n))
    ab[0, 1:] = -r/2
    ab[1, :]  = 1 + r
    ab[2, :-1] = -r/2
    for _ in range(Nt):
        ui  = u[1:-1]
        rhs = (1 - r) * ui
        rhs[1:]  += (r/2) * ui[:-1]
        rhs[:-1] += (r/2) * ui[1:]
        u[1:-1] = solve_banded((1, 1), ab, rhs)
    return x, u

波动方程

d’Alembert 与光锥

$$\boxed{\;u(x, t) = \frac{1}{2}\bigl[f(x - ct) + f(x + ct)\bigr] + \frac{1}{2c}\int_{x - ct}^{x + ct} g(s)\,ds.\;}$$

两条观察重写了"波"的含义：

$(x, t)$ 处的解只依赖依赖区间 $[x - ct, x + ct]$ 内的数据。区间外的扰动完全不可见——有限信号速度。
$f(x \pm ct)$ 是行波：刚体平移，分别向右、向左，速度都为 $c$。从静止释放的脉冲会分裂成两个等高脉冲，向相反方向以 $c$ 速度离开。

有限弦上的驻波

$$u(x, t) = \sum_{n=1}^\infty \bigl[a_n \cos(n\pi c\,t/L) + b_n \sin(n\pi c\,t/L)\bigr]\sin(n\pi x/L).$$

第 $n$ 个模态频率 $\omega_n = n\pi c/L$。$n=1$ 是基频；$n\geq 2$ 是泛音。这就是为什么单簧管和小提琴在同一音高下听感截然不同——泛音的相对幅度不同。

数值格式与 CFL 条件

$$u_j^{n+1} = 2u_j^n - u_j^{n-1} + \sigma^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n), \qquad \sigma = \frac{c\,\Delta t}{\Delta x}.$$

von Neumann 分析给出著名的 CFL 条件 $\sigma \leq 1$。几何解释非常优美：$(x, t)$ 处数值依赖区域——格式实际可能用到的网格点——必须包含物理依赖区域 $[x - ct, x + ct]$。一旦 $\sigma > 1$，物理信号比网格能传递的快，数值解必然爆炸。

波动方程：d’Alembert 分裂、驻波节点、时空光锥、有限差分快照、CFL 越界爆炸、特征曲线网。

上排：(a) 高斯脉冲分裂为两个半高脉冲反向离开；(b) 固定弦上 $n=3$ 的驻波——三段、两个内节点（红点处恒为零）；(c) 时空热图，虚线是光锥 $x = \pm ct$，外侧恒等于零。下排：(d) leapfrog 在 $\sigma=0.7$ 时干净地传播；(e) 同一格式 $\sigma=1.05$，半个时间单位内就出现数值湍流；(f) 特征曲线网 $x \pm ct = \text{常数}$——信息沿这些曲线传递。

Laplace 方程

与时间无关的平衡

$\nabla^2 u = u_{xx} + u_{yy} = 0$ 描述稳态分布——板内温度（边界温度长期保持后），引力势在质量外的形态，无旋流体的速度势。没有时间演化，只有边界数据。

两条结构性结果

$$\min_{\partial\Omega} u \leq u(\mathbf{x}) \leq \max_{\partial\Omega} u, \quad \forall \mathbf{x}\in\Omega.$$

内部不可能比边界最热点更热（也不可能比边界最冷点更冷）。两行证明的关键是均值性质：$u(\mathbf{x}_0)$ 等于以 $\mathbf{x}_0$ 为中心、完全包含于 $\Omega$ 的任意圆盘上 $u$ 的平均。如果内部某点取到最大值，则整个邻域内 $u$ 等于该最大值；连通性论证将其传播到边界。

唯一性：两个解差出一个边界为零的调和函数——由极值原理它必恒等于零。所以 Laplace BVP 恰有一个解。

数值方法：五点格式

$$u_{i, j} = \frac{1}{4}\bigl(u_{i+1, j} + u_{i-1, j} + u_{i, j+1} + u_{i, j-1}\bigr).$$

离散版的极值原理就是这个"取邻居平均"的语句——也直接给出算法：从任意初猜开始，反复用邻居平均替换内部值（Jacobi / Gauss-Seidel）。$N \times N$ 网格需要 $O(N^2)$ 次迭代；带最优松弛因子的 SOR 把它降到 $O(N\log N)$。实际工作中常用多重网格或稀疏直接解法，但迭代法是极值原理"暗示出来"的最自然算法。

Laplace 方程在单位正方形上的解：色图、等势线 + 热流方向、极值原理散点验证。

左：$\nabla^2 u = 0$ 的解，上边界为热高斯凸起，其他三边为零。热从上往下、向两侧渗透。中：等势线 $u = \text{常数}$ 与负梯度 $-\nabla u$（依 Fourier 定律为热流方向）。右：600 个内部随机采样点的 $u$ 值——全部位于边界最大值之下、最小值之上，极值原理得到直观验证。

Poisson 方程

加上源项：$\nabla^2 u = -f$。五点格式变为 $u_{i,j} = \tfrac{1}{4}(\text{邻居和}) + \tfrac{\Delta x^2}{4} f_{i,j}$。Green 函数把连续问题化为与基本解 $G(\mathbf{x}, \mathbf{y}) = -\frac{1}{2\pi}\ln|\mathbf{x} - \mathbf{y}|$（二维）的卷积——它是热方程高斯基本解的"源驱动型"对应物。

总结：方法地图

方程	类型	初/边值数据	数值方法	稳定性
$u_t = \alpha u_{xx}$	抛物型	$u(x, 0)$ + 2 个 BC	FTCS / Crank-Nicolson	$r \leq 1/2$（FTCS），无条件（CN）
$u_{tt} = c^2 u_{xx}$	双曲型	$u(x, 0),\ u_t(x, 0)$ + 2 个 BC	leapfrog	CFL $\sigma \leq 1$
$\nabla^2 u = 0$	椭圆型	仅 BC	Gauss-Seidel / SOR / 多重网格	始终稳定

面对一个新 PDE 的清单：

判型。线性二阶时直接算 $B^2 - AC$。
挑可行数据。双曲、抛物需要初值；椭圆不需要。
几何简单时选基（分离变量），复杂时离散化。
先看稳定性，再看精度。不稳定的格式根本谈不上精度。

PDE 的世界极其辽阔——非线性守恒律、流体力学、广义相对论、量子场论——我们才刚跨进门。但这三大正则方程是语法；其他都是词汇。

练习

概念题

为什么热方程"会忘掉初值"，而波动方程不会？请用 Fourier 模态语言说明。
证明：$\Delta = B^2 - AC$ 在 $(x, y)$ 平面任何旋转下不变。所以分类是几何性的而非坐标依赖。
用物理直觉解释为什么椭圆型问题不接受初值问题。

计算题

实现 Crank-Nicolson 求解热方程；通过同时减半 $\Delta t$ 与 $\Delta x$，验证二阶收敛。
数值探索波动方程在 $\sigma = 1$ 边界的行为：恰好等于 1 时格式有什么特殊？略大于 1 时呢？
直接代入验证 $G(x, t) = (4\pi\alpha t)^{-1/2}\,e^{-x^2/(4\alpha t)}$ 对 $t > 0$ 满足热方程。

编程题

在正方形上解 Laplace 方程：三边为 0，上边为 $\sin(\pi x/L)$。把迭代解与解析解 $u(x, y) = \sin(\pi x/L)\sinh(\pi y/L)/\sinh(\pi)$ 比较。
动画化三角脉冲初值的 d’Alembert 解；视觉验证它分裂为两个等高三角形。
重现 Gibbs 现象：用 $K = 5, 20, 100$ 个模态分别绘制方波的 Fourier 部分和，量化跳跃处的过冲幅度。

参考资料

Strauss, Partial Differential Equations: An Introduction, 2nd ed., Wiley (2007)
Haberman, Applied Partial Differential Equations, 5th ed., Pearson (2013)
Evans, Partial Differential Equations, 2nd ed., AMS (2010)
LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM (2007)
Courant, Friedrichs & Lewy, “On the partial difference equations of mathematical physics,” IBM J. Res. Develop. 11 (1967) —— CFL 条件原始论文
Press 等, Numerical Recipes, 3rd ed., Cambridge (2007), 第 19-20 章

系列导航


上一章	第十二章：边值问题
当前	第十三章：偏微分方程引论
下一章	第十四章：传染病模型与流行病学