3. 瞬态热传导显式与隐式时间推进对比

问题设定与离散框架

取一根长度 L = 0.1 m 的均质导热棒，热扩散率 α = 1.0×10⁻⁵ m²/s，初始温度均匀为 T₀ = 20 °C，两端固定为 T_L = T_R = 100 °C。控制方程为：

∂T/∂t = α ∂²T/∂x²

空间采用二阶中心差分，节点数 N = 21（含两端），Δx = L/(N−1) = 5 mm。时间推进分别使用显式欧拉与隐式欧拉，二者共享同一空间离散矩阵，仅时间项处理不同：

• 显式欧拉：Tⁿ⁺¹ = Tⁿ + Δt·α·(A·Tⁿ + b)，无需解线性方程组，但受稳定性约束。
• 隐式欧拉：(I − Δt·α·A)·Tⁿ⁺¹ = Tⁿ + Δt·α·b，每步需解三对角系统，但无条件稳定。

其中 A 为内部节点的二阶差分算子（除以 Δx²），b 承载边界温度贡献。显式格式的稳定性条件为 Δt ≤ Δx²/(2α) = 1.25 s，这是本节对比的核心参数。

显式与隐式推进的完整实现

下面给出可直接运行的 Python 实现，使用 NumPy 构建三对角矩阵并用 scipy.linalg.solve_banded 求解隐式步。代码同时输出两种格式在相同物理时间下的温度场与耗时，便于对比。

import numpy as np
from scipy.linalg import solve_banded
import time

# ---- 物理与网格参数 ----
L, alpha = 0.1, 1.0e-5
N = 21
dx = L / (N - 1)
T0, TL, TR = 20.0, 100.0, 100.0
t_end = 600.0            # 物理模拟总时长 (s)

# ---- 空间离散：内部节点 i=1..N-2 ----
n = N - 2
r = alpha / dx**2
main = -2.0 * r * np.ones(n)
off  = r * np.ones(n - 1)
A = np.diag(main) + np.diag(off, 1) + np.diag(off, -1)
b = np.zeros(n)
b[0]  += r * TL          # 左边界贡献
b[-1] += r * TR          # 右边界贡献

def init_field():
    T = np.full(N, T0)
    T[0], T[-1] = TL, TR
    return T

def explicit(dt):
    T = init_field()
    steps = int(round(t_end / dt))
    t0 = time.perf_counter()
    for _ in range(steps):
        Ti = T[1:-1]
        T[1:-1] = Ti + dt * (A @ Ti + b)
    return T, time.perf_counter() - t0, steps

def implicit(dt):
    T = init_field()
    steps = int(round(t_end / dt))
    # 构造带状矩阵 (I - dt*A)，上/下带宽各 1
    ab = np.zeros((3, n))
    ab[0, 1:]  = -dt * off      # 上对角线
    ab[1, :]   = 1.0 - dt * main
    ab[2, :-1] = -dt * off      # 下对角线
    t0 = time.perf_counter()
    for _ in range(steps):
        rhs = T[1:-1] + dt * b
        T[1:-1] = solve_banded((1, 1), ab, rhs)
    return T, time.perf_counter() - t0, steps

# ---- 对比运行 ----
dt_exp = 1.0                 # 满足 1.25 s 稳定上限
dt_imp = 10.0                # 隐式可取更大步长
Te, te, se = explicit(dt_exp)
Ti, ti, si = implicit(dt_imp)
print(f"显式: dt={dt_exp}s, 步数={se}, 耗时={te:.4f}s, 中点温度={Te[N//2]:.3f}°C")
print(f"隐式: dt={dt_imp}s, 步数={si}, 耗时={ti:.4f}s, 中点温度={Ti[N//2]:.3f}°C")

运行后应观察到：显式在 dt = 1.0 s 时稳定收敛，中点温度接近稳态值约 60 °C；若把 dt_exp 改为 1.5 s（超过 1.25 s 上限），温度会迅速振荡发散。隐式在 dt = 10 s 时仍稳定，但单步计算量更大。这正是两种格式在稳定性与效率上的本质差异。

稳定性、效率与精度的对比判读

下图把两种格式的关键差异组织为参数关系与结果判读视图：左侧是稳定性约束的来源，右侧是相同物理时长下的步数与耗时对比，底部给出精度判读依据。

判读要点：两种格式在稳态极限下给出相同温度场，差异集中在瞬态阶段与稳定性边界。显式格式的 Δt_max = 1.25 s 由网格与材料参数直接决定，网格加密时该上限按 Δx² 缩小，成本急剧上升；隐式格式虽无条件稳定，但 Δt 过大时一阶时间离散会引入明显的时间滞后误差，需通过步长收敛性检查确认精度可接受。