4. 显式与隐式时间推进格式

抛物型方程（以一维热传导方程 u_t = α u_xx 为典型代表）在空间离散之后，会得到一组常微分方程组 du/dt = A u + b，其中 A 来自空间差分算子。时间推进格式的选择，本质上就是决定“用哪一时刻的空间导数信息来更新下一时刻的解”。这一选择直接决定了格式的稳定性、时间步长上限以及每一步的计算代价，是抛物型方程数值求解中最核心的工程权衡之一。

为便于讨论，先约定符号：空间步长 Δx，时间步长 Δt，网格点 i = 0,1,…,N，时间层 n 表示 t_n = nΔt，u_iⁿ 表示第 n 时间层第 i 个网格点的数值解。空间二阶导数采用中心差分：(u_i+1 − 2u_i + u_i−1)/Δx²。定义网格比 r = αΔt/Δx²，它是贯穿本节所有稳定性判据的关键无量纲参数。

三种格式的构造与代数结构

前向欧拉（显式）：用第 n 层的空间导数推进到第 n+1 层，

u_i^{n+1} = u_i^n + r (u_{i+1}^n − 2u_i^n + u_{i−1}^n)

每个未知量 u_iⁿ⁺¹ 都可以由已知的第 n 层值直接算出，无需解线性方程组。稳定性分析（von Neumann 方法）给出放大因子 G = 1 − 4r sin²(θ/2)，要求 |G| ≤ 1，得到显式格式的稳定条件：

r = αΔt/Δx² ≤ 1/2  即  Δt ≤ Δx²/(2α)

注意 Δt 与 Δx² 成正比：空间网格加密一倍，时间步长必须缩小到四分之一。这是显式格式最致命的工程约束。

后向欧拉（隐式）：用第 n+1 层的空间导数，

u_i^{n+1} − r (u_{i+1}^{n+1} − 2u_i^{n+1} + u_{i−1}^{n+1}) = u_i^n

每一时间步需要求解一个三对角线性方程组。其放大因子 G = 1/(1 + 4r sin²(θ/2))，对任意 r > 0 都满足 |G| ≤ 1，因此是无条件稳定的。代价是每步一次三对角求解（可用 Thomas 算法，O(N) 复杂度），且时间精度只有一阶。

Crank-Nicolson（CN）：对第 n 层和第 n+1 层的空间导数取平均，

(u_i^{n+1} − u_i^n)/Δt = (α/2)[δ²u_i^{n+1} + δ²u_i^n]/Δx²

其中 δ² 表示中心二阶差分。整理后同样是一个三对角方程组，放大因子 G = (1 − 2r sin²(θ/2))/(1 + 2r sin²(θ/2))，对任意 r 无条件稳定，且时间精度为二阶。CN 在 r 较大时对高频分量衰减不足，初值或边界突变处可能出现非物理振荡，这是工程使用中需要留意的细节。

稳定性条件的工程含义

显式格式的 r ≤ 1/2 不是“建议值”，而是硬性约束。举例来说，若 α = 1.0×10⁻⁶ m²/s（典型金属热扩散系数量级），Δx = 1 mm，则 Δt ≤ (10⁻³)²/(2×10⁻⁶) = 0.5 s。若为了看清温度梯度把 Δx 加密到 0.1 mm，Δt 上限骤降到 0.005 s，时间步数增加 100 倍，而每步计算量也因网格点增加 10 倍而上升，总代价放大约 1000 倍。这正是显式格式在精细网格下“看起来简单、实际昂贵”的原因。

隐式格式则允许 Δt 远大于 Δx²/(2α)。但“无条件稳定”不等于“无条件准确”：时间步长过大时，一阶的后向欧拉会引入明显的数值耗散，把温度峰抹平；CN 虽然二阶，但在 r 很大时对初值高频分量的放大因子接近 −1，会产生衰减极慢的振荡。工程上的常见做法是：用隐式格式换取稳定性，同时把 Δt 控制在物理过程特征时间（如 αΔt/L² 不超过 0.1～1）以内，兼顾精度。

计算效率的实际权衡

不能只看“每步谁快”。显式每步只需一次向量化的差分更新，无需求解方程组，单步成本极低且极易并行；隐式与 CN 每步需要组装并求解三对角系统，单步成本约为显式的数倍到十余倍（取决于实现与问题规模）。但显式受 Δt ≤ Δx²/(2α) 限制，总步数随 Δx² 反比增长。因此存在一个交叉点：网格较粗、模拟时间较短时显式可能更划算；网格精细或需要长时间积分时，隐式/CN 的总耗时通常显著更低。

下面给出一个可直接用于估算的对比表（示例参数：α = 1×10⁻⁶ m²/s，L = 0.1 m，模拟总时长 T = 100 s，N 为空间网格点数）：

表中数值为按公式估算的示例，用于说明量级关系，并非某次实际运行的计时结果。实际选择时还应考虑边界条件的处理方式：显式格式的边界更新直接、隐式格式需把边界值并入右端项，CN 则要对边界做两层平均，若边界随时间剧烈变化，CN 的振荡风险更高，可改用后向欧拉或对前几步采用隐式启动。