4.2 热传导差分实验完整复盘

4.2 热传导差分实验完整复盘

三要素给了清单，本节把它变成一份完整实验记录：对象仍是第三章那根两端零温的杆。日志按五步走——用泰勒展开推显式格式并验相容阶、用冯·诺依曼分析推稳定域、在稳定域内验证收敛阶、故意违反稳定域看爆炸、最后请出无条件稳定的克兰克–尼科尔森格式收尾。每一步都有可复算的代码与预期输出，走完这份日志，4.1 的全部抽象概念都落到手感和数字上。

从泰勒展开推出格式

显式格式（学名 FTCS：时间向前、空间中心）的推导只用泰勒展开。网格 u_i^n 记 x_i = i·dx 处、t_n = n·dt 时刻的温度。时间导数用向前差商，空间导数用中心差商：

u_t ≈ (u_i^{n+1} − u_i^n)/dt，u_xx ≈ (u_{i+1}^n − 2u_i^n + u_{i−1}^n)/dx²。

代入 u_t = a²·u_xx，解出下一时刻：

u_i^{n+1} = u_i^n + r·(u_{i+1}^n − 2u_i^n + u_{i−1}^n)，其中 r = a²·dt/dx²。

格式的读法有物理味：某点的下一时刻温度 = 当前温度加上 r 倍的"与邻居的平均差"。温度比邻居热就降温，比邻居冷就升温——扩散的离散化身。相容阶由泰勒展开验收：把真解代入，时间差商的领先误差项含 dt/2（一阶），中心差商的领先误差项含 dx²/12（二阶），截断误差 O(dt) + O(dx²)。

冯·诺依曼分析：给稳定性定界

稳定性问题：舍入与初始扰动会不会被迭代放大？冯·诺依曼的思路是把扰动按傅里叶模式展开，追踪每个模式经一步迭代的放大因子。设扰动形式 gⁿ·e^(ikx_i)，代入格式：

g = 1 + r·(e^(ikdx) + e^(−ikdx) − 2) = 1 − 4r·sin²(k·dx/2)。

放大因子 g 是 r 与波数的函数。稳定性要求一切波数的 |g| 不超过 1。逐波数检查：g 在 sin² 取 1（最短波，k·dx = π）处最小，为 1 − 4r。要求 1 − 4r ≥ −1，即

r = a²·dt/dx² ≤ 1/2。

这条条件值得念出声：时间步长被空间步长平方锁死——网格加密一倍，步长必须缩到四分之一，总步数翻四倍，计算量按四次方暴涨。这是显式格式处理扩散问题的死穴，也是后面克兰克–尼科尔森登场的理由。注意条件只卡 r 与波长的关系，与杆长、初值大小无关——纯格式性质，纯泰勒展开加复指数代数，纸笔十分钟。

实验一：稳定域内验证收敛阶

固定 r = 0.4（安全），三套网格对表（4.1 已写过框架，这里补完整跑法与读数）。以 3.1 的级数解为标尺：

import numpy as np

L, a, T_end, r = 1.0, 1.0, 0.1, 0.4
series = lambda x, t: sum(8/(n**3*np.pi**3)*np.sin(n*np.pi*x/L)
                          * np.exp(-a**2*(n*np.pi/L)**2*t) for n in range(1, 4000, 2))

def run(nx):
    x = np.linspace(0, L, nx); dx = x[1]-x[0]; dt = r*dx**2/a**2
    u = x*(L-x)
    for _ in range(int(round(T_end/dt))):
        u[1:-1] += r*(u[2:] - 2*u[1:-1] + u[:-2])
    return np.abs(u - series(x, T_end)).max()

errs = [run(nx) for nx in [11, 21, 41, 81]]
for i, e in enumerate(errs):
    rate = "" if i == 0 else f"  阶数约 {np.log2(errs[i-1]/e):.2f}"
    print(f"nx={11*2**i:3d}  误差 {e:.3e}{rate}")
​

预期输出：相邻误差比稳定在 4 附近，log₂ 后阶数约 2——空间二阶实测兑现（时间误差被 r 联动控制，同阶）。阶数一旦只有 1，先查边界格式；一旦乱跳，先查时间终点与步数是否整除对齐——两大调试经验顺手记下。

实验二：越界看爆炸

把 r 提到 0.6，其余不动。预期现象分阶段：前几十步解还能看，随后高频锯齿（最短波长的棋盘模式）从舍入噪声里长出来，指数增厚，几十步后数值溢出。这正是放大因子预言的：r 大于 1/2 时 g 在最短波处低于 −1，棋盘模式每步反号且幅值增长。理论与现象的对账精确到"哪种波形先炸"——不是随便乱炸，是 k·dx = π 的模式先炸，因为它的 g 最负。想亲眼看模式选择，打印爆炸前的解减去均值，会看到规整的交替正负棋盘。

r_bad = 0.6
x = np.linspace(0, L, 21); dx = x[1]-x[0]; dt = r_bad*dx**2/a**2
u = x*(L-x)
for step in range(1, 200):
    u[1:-1] += r_bad*(u[2:] - 2*u[1:-1] + u[:-2])
    if not np.all(np.isfinite(u)):
        print(f"第 {step} 步溢出  dt = {dt:.2e}")
        break
​

实验三：克兰克–尼科尔森修复死穴

把时间导数改在 n 与 n+1 两层的平均位置取值，空间项也取两层平均：

u_i^{n+1} − u_i^n = (r/2)·[(u_{i+1} − 2u_i + u_{i−1})^{n+1} + (u_{i+1} − 2u_i + u_{i−1})^n]。

代价是每步解一个三对角线性系统（追赶法 O(N)，廉价）；回报是放大因子变为 (1 − 2r·sin²)/(1 + 2r·sin²)，恒小于 1——无条件稳定，步长不再被网格平方锁死，且时间精度升至二阶。数值上把它与显式格式同台对比：同样网格精度下，克兰克–尼科尔森可以用大得多的步长。一个值得记录的残余问题：大步长下高频分量的放大因子趋于 −1，初始棱角的高频尾巴衰减缓慢，解会带轻微的振荡拖尾——修正手段是把时间导数略向隐式偏置（所谓 Rannacher 平衡步）。这个细节在期权定价（第六章 Black-Scholes 数值解）里是必踩的坑，那里先记一笔。

⚠️ 常见坑：无条件稳定不等于无条件精确。步长放到很大时稳定性无恙，时间误差照常增长——精度的取舍永远要重新对表，"稳定"只是入场券。

本节要点

• FTCS 格式：泰勒展开两条差商拼成，物理读法是"向邻居均值回归"，截断 O(dt) 加 O(dx²)；
• 冯·诺依曼分析：追踪每个傅里叶模式的放大因子，r 不超过二分之一，纸笔十分钟定生死；
• 平方死穴：扩散显式格式的计算量随加密按四次方涨，无条件稳定的隐式格式因此成为工程标配；
• 爆炸有模式：违反条件时最短波的棋盘模式最先失控，现象与理论精确对账；
• 克兰克–尼科尔森：两层平均换来无条件稳定与二阶时间精度，代价是每步解三对角系统；
• 实验纪律：阶数对表、越界复爆、大步长查拖尾——三组对照是一切格式验收的模板。

差分把导数换成差商，假设是网格规整。真实结构的几何从不规整——下一节的有限元换一条路线：在任意三角形拼出的区域上解弱形式，复杂几何从此不再是障碍。