4. 显式格式的编程实现与验证

本节把前面推导的FTCS（Forward-Time Central-Space，时间前向、空间中心）离散格式真正落到代码上，并用一个具有解析解的金属棒瞬态导热案例来检验实现是否正确。学习目标不是"跑通一段程序"，而是建立一条可复现的验证链：明确物理问题与边界条件 → 按FTCS公式组装更新循环 → 用解析解逐点比对 → 用误差范数给出定量结论 → 用时空演化图观察物理合理性。只有这五步都成立，才能说求解器"可信"。

一、物理问题与解析解基准

取一根长度 L = 1 m 的均质金属棒，热扩散系数 α = 1.0×10⁻⁴ m²/s（量级接近不锈钢），初始温度分布取正弦模态：

T(x,0) = T₀ + A·sin(πx/L)，T₀ = 20 ℃，A = 100 ℃

两端保持恒定温度 T(0,t) = T(L,t) = T₀（第一类边界条件）。该问题的解析解为：

T(x,t) = T₀ + A·sin(πx/L)·exp(−α π² t / L²)

选择正弦初值而非阶跃初值，是因为正弦模态是热方程的特征函数，解析解光滑、无奇点，最适合做代码验证；阶跃初值在 x=0 处导数不连续，会引入额外的离散误差，不适合作为第一轮验证对象。注意温度以 T₀ 为基准，程序中可令 u = T − T₀ 简化计算，最后再加回 T₀。

二、离散参数与稳定性约束

FTCS 的更新公式为：

u[i]^(n+1) = u[i]^n + r·(u[i+1]^n − 2u[i]^n + u[i−1]^n)，r = α·Δt/Δx²

显式格式的稳定性条件是 r ≤ 0.5。取 Nx = 51 个节点（Δx = 0.02 m），若取 r = 0.4，则 Δt = r·Δx²/α = 0.4×4×10⁻⁴/1×10⁻⁴ = 1.6 s。总模拟时长取 t_end = 200 s，约需 125 步。这个参数组合既满足稳定性，又让时间推进步数适中，便于观察。若误取 r = 0.6，数值解会出现非物理的振荡并迅速发散——这正是验证环节要捕捉的典型错误。

三、Python 实现要点

核心循环用 NumPy 向量化切片实现，避免逐点 Python 循环。关键写法是：u[1:-1] = u[1:-1] + r*(u[2:] - 2*u[1:-1] + u[:-2])，其中 u[2:] 对应 u[i+1]，u[:-2] 对应 u[i−1]。边界节点 u[0] 与 u[-1] 不参与更新，保持为 0（即 T₀）。

import numpy as np
import matplotlib.pyplot as plt

L, alpha, T0, A = 1.0, 1e-4, 20.0, 100.0
Nx = 51
dx = L / (Nx - 1)
x = np.linspace(0, L, Nx)
r = 0.4
dt = r * dx**2 / alpha
t_end = 200.0
nsteps = int(round(t_end / dt))

u = A * np.sin(np.pi * x / L)   # 初始条件（相对温度）
u[0] = u[-1] = 0.0              # 边界条件

snapshots = {0: u.copy()}
for n in range(1, nsteps + 1):
    u[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
    t = n * dt
    if abs(t - 50) < dt/2 or abs(t - 100) < dt/2 or abs(t - 200) < dt/2:
        snapshots[round(t)] = u.copy()

def exact(x, t):
    return A * np.sin(np.pi * x / L) * np.exp(-alpha * np.pi**2 * t / L**2)

for t, un in snapshots.items():
    err = un - exact(x, t)
    Linf = np.max(np.abs(err))
    L2 = np.sqrt(dx * np.sum(err**2))
    print(f"t={t:3d}s  Linf={Linf:.4e}  L2={L2:.4e}")

代码中有两个容易出错的细节。第一，nsteps 用 round 而非 int 截断，避免因浮点误差少走一步；第二，快照判断用 abs(t - 目标) < dt/2，因为 n*dt 未必精确等于整数秒。若直接用 t == 50 判断，可能永远不触发。

四、误差范数与收敛性判读

在 t = 50、100、200 s 三个时刻分别计算 L∞ 与 L2 误差。以本节参数为例，合理预期是：L∞ 误差在 10⁻² ℃ 量级，且随时间增大而略有增长（因为时间累积误差），但不会发散。若出现 L∞ 达到 10¹ 以上或随时间爆炸，几乎可以断定 r > 0.5 或边界处理写错。

更严格的验证是网格收敛性检验：固定 r = 0.4，将 Nx 从 51 加密到 101、201，观察 L∞ 是否按约 4 倍（Δx² 主导）下降。由于 FTCS 的时间误差是一阶、空间误差是二阶，在 r 固定时 Δt ∝ Δx²，时间误差也退化为二阶，因此整体应呈现接近二阶收敛。若加密后误差不降反升，说明实现中存在索引错位或边界泄漏。

五、时空演化图的观察要点

用 plt.imshow 或 plt.pcolormesh 绘制 T(x,t) 的二维色图，横轴 x、纵轴 t、色标为温度。正确的图应呈现：初始时刻中间高、两端低的半正弦拱形；随时间推移，峰值沿时间轴指数衰减并趋于 T₀；两端始终为 T₀ 的竖直等温带；整个场无棋盘状振荡、无负温度。若出现沿 x 方向的高频锯齿，通常是 r 接近 0.5 时的临界振荡；若出现负值，则已越过稳定性边界。