4. 代码验证

时间步收敛性研究的设计

时间步收敛性分析的核心是：固定空间离散（网格、空间格式、边界条件不变），只改变时间步长 Δt，观察某个标量输出量 Q（如某时刻探针温度、总能量、峰值位移）如何随 Δt 减小而变化。若时间积分格式为 p 阶精度，则误差 E(Δt) = |Q(Δt) − Q_ref| 应满足 E ≈ C·Δt^p，在双对数坐标下呈斜率 −p 的直线。

参考解 Q_ref 的获取有两种稳妥做法：一是用极小的 Δt（如比最小测试步长再小两个数量级）计算，二是用两种不同阶格式交叉验证。本例采用第一种，并明确标注 Q_ref 为数值参考解而非解析解。

稳定性检查：不只是"算得动"

显式格式受 CFL 条件约束。以一维热传导方程为例，显式 FTCS 格式的稳定性要求 r = α·Δt/Δx² ≤ 0.5。当 r 超过 0.5 时，数值解会出现非物理振荡并迅速发散。因此时间步收敛性研究必须同时覆盖稳定区和不稳定区，才能区分"精度不足"与"稳定性失效"两种截然不同的误差来源。

稳定性检查的实操方法：在收敛性扫描中额外加入一个 r = 0.6 的算例，观察解是否出现振荡放大。若该算例发散而 r ≤ 0.5 的算例正常收敛，则说明稳定性边界判断与理论一致，代码的稳定性行为正确。

收敛性扫描的完整流程

下面给出一个可直接运行的时间步收敛性分析脚本。它求解一维瞬态热传导方程，使用显式 FTCS 格式，对一系列 Δt 计算终点时刻的温度场，并以最小步长解为参考计算误差和观测收敛阶。

import numpy as np

# ---- 问题设置 ----
L = 1.0            # 长度 [m]
alpha = 1e-4       # 热扩散率 [m^2/s]
nx = 51            # 空间节点数（固定）
dx = L / (nx - 1)
x = np.linspace(0, L, nx)

T_left, T_right = 100.0, 0.0   # 边界条件
T_init = 0.0                    # 初始条件
t_end = 50.0                    # 终点时刻 [s]

def solve_ftcs(dt):
    """显式FTCS求解一维热传导，返回终点温度场"""
    nt = int(round(t_end / dt))
    dt_actual = t_end / nt          # 调整使恰好到达t_end
    r = alpha * dt_actual / dx**2
    if r > 0.5:
        return None, r, dt_actual   # 超出稳定域，标记
    T = np.full(nx, T_init)
    T[0], T[-1] = T_left, T_right
    for _ in range(nt):
        T_new = T.copy()
        T_new[1:-1] = T[1:-1] + r * (T[2:] - 2*T[1:-1] + T[:-2])
        T = T_new
    return T, r, dt_actual

# ---- 时间步序列 ----
dt_list = [0.5, 0.25, 0.125, 0.0625, 0.03125]

# 参考解：极小步长
T_ref, r_ref, dt_ref = solve_ftcs(0.005)
print(f"参考解: dt={dt_ref:.5f}s, r={r_ref:.4f}")

# ---- 扫描与误差计算 ----
print(f"\n{'dt':>8} {'r':>8} {'E_max':>12} {'p_obs':>8}")
results = []
for dt in dt_list:
    T, r, dt_a = solve_ftcs(dt)
    if T is None:
        print(f"{dt:>8.5f} {r:>8.4f}   —— 超出稳定域(r>0.5)")
        continue
    E = np.max(np.abs(T - T_ref))
    results.append((dt_a, E, r))
    if len(results) >= 2:
        dt1, E1, _ = results[-2]
        dt2, E2, _ = results[-1]
        p_obs = np.log(E1/E2) / np.log(dt1/dt2)
        print(f"{dt_a:>8.5f} {r:>8.4f} {E:>12.3e} {p_obs:>8.3f}")
    else:
        print(f"{dt_a:>8.5f} {r:>8.4f} {E:>12.3e} {'—':>8}")

# ---- 稳定性边界验证 ----
print("\n--- 稳定性边界检查 ---")
for dt_test in [0.4, 0.5, 0.6]:
    T, r, _ = solve_ftcs(dt_test)
    status = "稳定" if T is not None else "超出稳定域"
    print(f"dt={dt_test:.2f}s, r={r:.4f} → {status}")

运行后预期看到：随着 Δt 减半，E_max 约减半（p_obs ≈ 1），因为显式 FTCS 的时间精度是一阶。若 p_obs 稳定在 1.0 附近，说明时间积分实现正确；若 p_obs 偏离或随步长变化，需检查边界处理、源项离散或迭代收敛容差。

结果判读与常见陷阱

时间步收敛性分析中最容易踩的坑有三类：

• 空间误差未排除：若网格本身未收敛，减小 Δt 时误差可能被空间误差"地板"截断，p_obs 趋近于 0。应先用网格收敛性分析确认空间离散已足够细。
• 非线性迭代容差过松：隐式格式中若每步迭代未收敛到足够精度，时间精度会被迭代误差掩盖，表现为 p_obs 低于理论阶。
• 稳定性边界误判：r 恰好接近 0.5 时，解可能暂时稳定但误差放大，导致收敛阶异常。建议测试步长覆盖 r = 0.1 ~ 0.45 的范围。

实操检查清单