有限差分与热传导：用 Python 搭一个可验证的数值仿真

从一个物理问题开始#

考虑一块正方形薄板，四边温度固定为 0，初始时刻中心区域温度较高。我们想知道温度如何随时间衰减。

二维热传导方程可以写成：

∂T/∂t = α (∂²T/∂x² + ∂²T/∂y²)

其中 T 是温度，α 是热扩散率。这个方程没有一般初边值条件下的解析解，但它的正弦模态有精确解，正好可以用来验证代码。

显式差分：把导数换成差值#

在均匀网格上，用中心差分近似二阶空间导数，用前向差分近似时间导数：

T[i,j]^(n+1) = T[i,j]^n + r * (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1] - 4T[i,j])

其中：

r = α * Δt / Δx²

这个格式叫 FTCS（Forward Time, Central Space）。它实现简单，但不是无条件稳定的。对二维正方形网格，稳定条件为：

r ≤ 1/4

如果时间步取得太大，数值解会振荡并发散。这个限制不是实现细节，而是格式本身的性质。

稳定性与精度是两件事#

稳定只代表“不会炸”，不代表“算得准”。工程报告里只写“没有报错”是不够的。

完整实现#

下面的代码用 NumPy 的切片做向量化更新，避免逐点循环。边界始终保持为 0。

import numpy as np
 
# 物理与网格参数
L = 1.0
alpha = 1e-4
nx = ny = 81
dx = L / (nx - 1)
 
# 稳定性：二维显式格式要求 r <= 0.25
r = 0.2
dt = r * dx**2 / alpha
t_end = 1000.0
steps = int(t_end / dt)
 
x = np.linspace(0.0, L, nx)
y = np.linspace(0.0, L, ny)
X, Y = np.meshgrid(x, y, indexing="ij")
 
# 初始条件：满足零边界的一阶正弦模态
T = np.sin(np.pi * X) * np.sin(np.pi * Y)
T[0, :] = T[-1, :] = 0.0
T[:, 0] = T[:, -1] = 0.0
 
for _ in range(steps):
    T_new = T.copy()
    T_new[1:-1, 1:-1] = T[1:-1, 1:-1] + r * (
        T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2]
        - 4.0 * T[1:-1, 1:-1]
    )
    T = T_new
 
print("steps:", steps, "dt:", dt)
print("center temperature:", T[nx // 2, ny // 2])

这段代码有三个值得注意的点：

• T_new = T.copy() 保留了边界值，避免边界被内部更新覆盖。
• 更新只作用于 [1:-1, 1:-1]，边界永远不参与计算。
• r = 0.2 留出了安全余量，方便后续改参数。

用解析解验证#

初始条件取 sin(πx)sin(πy) 时，精确解为：

T(x,y,t) = sin(πx) sin(πy) exp(-2 α π² t)

在 t = 1000 时，可以计算数值解与解析解的相对 L2 误差：

exact = np.sin(np.pi * X) * np.sin(np.pi * Y) * np.exp(-2 * alpha * np.pi**2 * t_end)
error = np.sqrt(np.mean((T - exact) ** 2)) / np.sqrt(np.mean(exact**2))
print("relative L2 error:", error)

在 81×81 网格、r = 0.2 的设置下，误差通常在百分之一量级。若误差明显偏大，优先检查：

1. 时间步数是否取整导致实际结束时间偏离 t_end。
2. 边界是否在迭代中被意外修改。
3. meshgrid 的 indexing 与数组索引是否一致。

网格无关性：别只看一张图#

把 nx 从 41 加密到 81、161，保持 r 不变，观察中心温度或 L2 误差。二维二阶格式在网格加密一倍时，误差理论上是四分之一。实际曲线会在某个尺度后变平，原因可能是：

• 时间误差开始主导；
• 边界附近的离散误差累积；
• 浮点误差接近机器精度。

找到“误差不再明显下降”的区间，比盲目追求最细网格更有意义。

工程实现的几个提醒#

• 单位统一：α 的单位是 m²/s，长度用米，时间用秒。混用毫米会让 r 差出六个数量级。
• 先算稳定性再跑：把 dt 写成 r * dx**2 / alpha，而不是手填一个看起来合理的数。
• 保存配置：网格数、r、alpha、结束时间、随机种子，全部写进结果文件。
• 不要用渲染图代替误差：温度云图好看，不代表数值正确。

可以继续做的扩展#

1. 把显式格式换成 Crank–Nicolson，比较同样精度下的时间步需求。
2. 加入非均匀材料，让 α 随空间变化。
3. 把边界条件换成对流换热，观察边界层。
4. 用稀疏矩阵组装，交给 scipy.sparse.linalg 求解隐式格式。

先把这个小闭环跑通，再扩展到真实几何，路线会稳得多。