6.2 案例研究：模拟与验证

6.2 案例研究：模拟与验证

本节摘要：以一维热传导方程为完整案例，演示三道验证关卡：守恒量检查（总热量是否闭合）、网格收敛研究（实测收敛阶是否与理论一致）、解析基准对照（与误差函数精确解比对）。末尾给出数值结果报告的标准模板。

一、案件设定

一维热传导：杆长 1，初值为中心高温脉冲，两端绝热。用隐式欧拉（4.3 节：扩散问题显式步长受 h 平方枷锁，隐式免于此难）离散：

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

def solve_heat(N=200, dt=1e-3, T_final=0.1, diffusion=0.01):
    h = 1.0 / (N - 1)
    x = np.linspace(0, 1, N)
    u = np.exp(-200 * (x - 0.5)**2)          # 初始温度脉冲
    r = diffusion * dt / h**2
    total_heat0 = u.sum() * h                 # 初始总热量（梯形近似）
    # 绝热边界：内部点离散，两端 du/dx = 0
    n = N - 2
    A = diags([np.full(n, -r), np.full(n, 1 + 2*r), np.full(n, -r)],
              [-1, 0, 1], format='csc')
    steps = int(T_final / dt)
    for _ in range(steps):
        u_int = spsolve(A, u[1:-1])           # 隐式步：解线性系统（第二章工具）
        u = np.concatenate([[u[0]], u_int, [u[-1]]])
    return x, u, total_heat0, u.sum() * h
​

程序跑通，图形上脉冲平滑展开——一切"看起来"正常。现在开始三关验证。

二、第一关：守恒量检查

绝热边界下总热量必须守恒。这是免费的正确性证词：不需要知道真解，只需检验模拟自身的物理约束：

import numpy as np

x, u, heat0, heat_final = solve_heat(N=201, dt=1e-3, T_final=0.1)
print(f"初始总热量: {heat0:.12f}")
print(f"末态总热量: {heat_final:.12f}")
print(f"相对漂移:   {abs(heat_final-heat0)/heat0:.2e}")
# 隐式格式对线性守恒律天然保持，漂移应在 1e-12 量级（纯舍入级）
​

若漂移达到 1e-3 量级，几乎必然是边界处理或装配出 bug——守恒检查是最灵敏的实现错误探测器。更多守恒证词：动量、能量（4.2 节的谐振子体检）、质量（反应系统组分和）、概率归一化（MCMC 的建议分布对称性）。任何 PDE/ODE 仿真上手第一件事：找到系统的守恒量并把它打印出来。

三、第二关：网格收敛研究

把网格加密一倍，结果变化多少？这个问题的答案决定结果的可信位数：

import numpy as np

def center_value(N, dt):
    x, u, _, _ = solve_heat(N=N, dt=dt)
    return u[N // 2]                          # 监测点：杆中心温度

h_prev = None
for N, dt in [(101, 2e-3), (201, 1e-3), (401, 5e-4), (801, 2.5e-4)]:
    v = center_value(N, dt)
    h = 1.0 / (N - 1)
    if h_prev is not None:
        ratio = (v_prev - v)  # 相邻两层解的差
        print(f"N={N:4d}  h={h:.5f}  中心温度 {v:.8f}  与上层差 {abs(ratio):.2e}")
    v_prev, h_prev = v, h
​

判读标准：相邻差按理论阶衰减（空间二阶则差值约缩 4 倍），且差值本身小于交付精度要求。实测阶数对不上理论阶，等价于实现有 bug——这是比任何单元测试都强的正确性证明，因为它同时校验了离散、装配、求解器全链路。若差值迟迟不降（收敛到错误的极限），说明离散公式本身错了。

图 6.2-1 网格收敛研究判读图（示意）

四、第三关：解析基准对照

热传导有误差函数精确解，初始脉冲演化可解析表达。有精确解的场景（通常通过特殊初值构造）是金标准：

import numpy as np
from scipy.special import erf

def analytic(x, t, diffusion=0.01, a=200.0):
    # 高斯脉冲在无限长杆上的解析解（卷积扩散核）
    s2 = a / (a + 4 * diffusion * t * a * a / 1)  # 简化：直接用扩散后的方差
    var = 1.0/(2*a) + 2 * diffusion * t
    return 1.0/np.sqrt(2*np.pi*var) * np.exp(-(x-0.5)**2/(2*var))

x, u_num, _, _ = solve_heat(N=801, dt=2.5e-4, T_final=0.1)
u_exact = analytic(x, 0.1)
err = np.abs(u_num - u_exact).max()
print(f"与解析解的最大偏差: {err:.3e}")
# 偏差应与理论截断误差同量级；若大几个数量级，回到第二关查收敛阶
​

常用基准库思路（纯文字指引，按关键词可查）：热传导误差函数解、波动方程 d'Alembert 解、泊松方程多项式特解（MS 法构造）、Burgers 方程 Cole-Hopf 变换解。没有解析解时，用高精度参考解替代——把网格和时间步加密到极限算一次当真解（本例的 801 网格结果就可以当粗网格的参考）。

五、结案模板：数值结果的标准报告

三关全过后，交付物按此模板撰写（每项都对应前文的某次取证）：

1. 数值结论：中心温度 0.1123（不是"约 0.11"这种含糊话）
2. 误差来源清单：空间离散二阶截断、时间一阶、浮点舍入（通常可忽略）、模型参数不确定（如扩散系数 ±5%）
3. 精度证据：网格收敛研究显示相邻层差 3e-4 且按理论阶衰减；守恒漂移 2e-12；与解析基准最大偏差 4e-4
4. 可信度陈述：数值误差贡献约 3e-4，参数不确定性贡献约 5e-3——结论的瓶颈在参数测量而非计算
5. 复现信息：随机种子（如有）、软件版本、关键参数

第 4 条是全模板的点睛之笔：它把"计算误差"与"模型/参数误差"放在同一杆秤上，明确告诉决策者改进哪里才有收益——这正是 5.2 节敏感性分析精神的落点。

⚠️ 常见坑：只做"图形看起来合理"的目视检查。对称、光滑、量级正常只能排除低级错误，收敛阶错误、边界条件用错等深层 bug 在图形上完全不显形。

本节要点回顾

• 守恒检查免费且灵敏：不知道真解也能做，实现错误的第一探测器
• 收敛阶是强正确性证明：实测阶对上理论阶，等于校验了全链路
• 收敛不等于正确：可能收敛到错误极限，第三关的基准对照不可省
• 三关递进：守恒（快）→ 收敛研究（中）→ 基准对照（慢但金标准）
• 报告第 4 条：把计算误差与参数误差同秤比较，指出瓶颈所在

下一节把前五章的所有案件浓缩成一本症状驱动的排错速查手册。