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

一旦某个量依赖于多个变量，常微分方程（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 章的边值问题，以及多元微积分基础（梯度、Laplacian 算子、链式法则）。

从 ODE 到 PDE#

在 ODE 中，未知函数 $$y(t)$$ 仅依赖于单个变量。每个解由有限个常数参数化——数量等于微分方程的阶数——而一个初始条件即可从中选出唯一解。

PDE 则要求解一个多变量函数。要确定它，必须在整个低维曲面（二维情形下是一条曲线，三维则是曲面）上给出数据，其解空间本质上是无限维的。同一个方程，对某种类型的数据可能是适定的，但换一种数据却可能不适定——这种现象在 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{常数}).

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

\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$$ 的半衰期与 $$1/n^2$$ 成正比。初始条件中的尖锐特征几乎瞬间被抹平，而平滑特征则缓慢消退。
边界条件决定了基函数的选择。Dirichlet 条件对应正弦函数，Neumann 对应余弦，周期性边界则对应复指数函数。这些本征函数彼此正交，因此系数可通过内积计算。这正是 Laplacian 算子的谱分解。

左图： $$u(x, t)$$ 在六个时刻的快照——一个分段常数的初始轮廓在不到一秒内被平滑，并逐渐趋近于平衡态 $u \equiv 0$ 。中图：同一解的时空热力图——明亮结构随时间向上迅速褪色。右图：Fourier 模态的对数尺度衰减曲线；模态 $$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^{ik j \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$$ 向右和向左传播。一个静止释放的脉冲会分裂为两个振幅减半的子脉冲，彼此远离。

有限弦上的驻波#

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.

这意味着区域内部不可能比边界最热点更热，也不可能比最冷点更冷。其证明仅需两行，核心是均值性质： $u(\mathbf{x}_0)$ 等于以 $\mathbf{x}_0$ 为中心、完全包含于 $\Omega$ 的任意圆盘上 $$u$$ 的平均值。若存在内部最大值，则它必须等于邻域内所有点的值，再通过连通性可推出该最大值必延伸至边界。

唯一性：若两个解具有相同的边界数据，则它们的差在边界上为零。由极值原理可知，该差函数在整个区域内恒为零。因此，Laplace 边值问题恰有唯一解。

数值方法：五点差分格式#

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$ （根据傅里叶定律，即热流方向）叠加显示。右图：在区域内随机选取 600 个点，绘制其 $$u$$ 值——所有点均严格位于边界最大值之下、最小值之上，直观验证了极值原理。

Poisson 方程#

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 个边界条件	FTCS / Crank-Nicolson	$r \leq 1/2$ （FTCS），无条件稳定（CN）
$u_{tt} = c^2 u_{xx}$	双曲型	$u(x, 0),\ u_t(x, 0)$ + 2 个边界条件	leapfrog	CFL 条件 $\sigma \leq 1$
$\nabla^2 u = 0$	椭圆型	仅需边界条件（无初值）	Gauss-Seidel / SOR / 多重网格	始终稳定

面对一个新的 PDE，可按以下清单逐步处理：

判断类型：若是线性二阶方程，直接计算 $$B^2 - AC$$ 。
选择合适的数据：双曲型和抛物型问题需要初始条件；椭圆型则不需要。
选择求解策略：几何简单时用分离变量法；复杂几何时采用离散化方法。
优先检查稳定性：在追求精度之前，务必确保数值格式稳定——不稳定的格式，精度毫无意义。

PDE 的世界浩瀚无边——非线性守恒律、流体力学、广义相对论、量子场论……我们才刚刚跨过门槛。但这三个经典方程构成了该领域的语法，其余一切不过是词汇的扩展。

练习题#

概念题

为何热方程会“遗忘”初始条件，而波动方程不会？请用 Fourier 模态的语言解释。
证明在 $$(x, y)$$ 平面的任意旋转下，判别式 $\Delta = B^2 - AC$ 保持不变，从而说明 PDE 的分类是几何性质，与坐标系选择无关。
从物理角度解释，为何椭圆型方程无法提出初值问题。

计算题

实现 Crank-Nicolson 方法求解热方程，并通过同时将 $\Delta t$ 和 $\Delta x$ 减半来验证其二阶收敛性。
对波动方程，数值探究 $\sigma = 1$ 附近的边界行为：当 $\sigma$ 恰好等于 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 部分和，并测量跳跃点处的过冲幅度。

下一步#

下一章把 ODE 应用到一个具体而紧迫的领域——传染病动力学。SIR 模型、再生数 $$R_0$$ 、群体免疫阈值，这些 2020 年家喻户晓的概念都建立在几个简单 ODE 之上。下一章不仅讲模型，更展示如何从数据估参、为什么相同模型在不同地区表现差异巨大。

参考文献#

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 章