3. 显式向前欧拉格式推导与稳定性分析

一维热方程的标准形式为

∂u/∂t = α ∂²u/∂x²，  0 < x < L，t > 0

其中 α > 0 为热扩散系数（单位 m²/s），u(x,t) 为温度。本节从最直接的差分离散出发，推导 FTCS（Forward-Time Central-Space，时间向前、空间中心）显式格式，并用冯·诺依曼方法给出稳定性条件，最后通过一个可手算的算例观察稳定与不稳定两种现象。

一、网格与差分算子

取空间步长 Δx，节点 x_i = iΔx（i = 0,1,…,N），时间步长 Δt，时刻 t_n = nΔt。记数值解 U_iⁿ ≈ u(x_i, t_n)。时间方向用一阶向前差商，空间方向用二阶中心差商：

∂u/∂t ≈ (U_iⁿ⁺¹ − U_iⁿ) / Δt

∂²u/∂x² ≈ (U_i+1ⁿ − 2U_iⁿ + U_i−1ⁿ) / Δx²

代入热方程并整理，得到 FTCS 格式：

U_iⁿ⁺¹ = U_iⁿ + r (U_i+1ⁿ − 2U_iⁿ + U_i−1ⁿ)，  r = αΔt / Δx²

这里 r 称为网格傅里叶数（也称扩散数）。该格式的截断误差为 O(Δt) + O(Δx²)，时间一阶、空间二阶。它的最大优点是每一步只需已知上一层的三个相邻值，无需解线性方程组，实现极为简单；代价是 Δt 受严格限制。

二、冯·诺依曼稳定性分析

把误差也看作满足同一线性差分方程的量，设误差模态为 ε_iⁿ = Gⁿ e^{jβiΔx}，其中 j 为虚数单位，β 为波数，G 为放大因子。代入 FTCS 格式并约去公因子 Gⁿe^{jβiΔx}：

G = 1 + r (e^{jβΔx} − 2 + e^{−jβΔx}) = 1 − 2r(1 − cos βΔx) = 1 − 4r sin²(βΔx/2)

稳定性要求对所有波数都有 |G| ≤ 1。由于 sin²(βΔx/2) ∈ [0,1]，G 的取值范围是 [1 − 4r, 1]。要保证 G ≥ −1，必须 1 − 4r ≥ −1，即

r = αΔt / Δx² ≤ 1/2，  等价于 Δt ≤ Δx² / (2α)

这就是 FTCS 格式的稳定性条件。它说明：空间网格越细（Δx 越小），允许的时间步长按平方缩小，计算量急剧上升。当 r > 1/2 时，高频误差模态的 |G| > 1，误差会逐步放大，数值解出现非物理的振荡并最终溢出。

三、算例：同一初值下的两种时间步长

取 α = 1 m²/s，杆长 L = 1 m，空间步长 Δx = 0.1 m（N = 10 个区间），初值为中间高、两端低的三角温度分布，两端边界固定为 0 ℃。由稳定性条件 Δt ≤ Δx²/(2α) = 0.01/2 = 0.005 s。下面分别取 Δt = 0.004 s（r = 0.4，稳定）与 Δt = 0.006 s（r = 0.6，不稳定）做对比。

稳定情形下，每步更新为 U_iⁿ⁺¹ = 0.4U_i−1ⁿ + 0.2U_iⁿ + 0.4U_i+1ⁿ，系数全为正且和为 1，相当于对邻点做加权平均，温度峰值单调下降、分布逐渐变平，符合物理扩散规律。不稳定情形下系数变为 0.6U_i−1ⁿ − 0.2U_iⁿ + 0.6U_i+1ⁿ，中心点系数为负，相邻节点出现正负交替的锯齿振荡，振幅随步数指数增长，几步之后数值即失去物理意义。

一个便于课堂演示的判据是：显式格式的更新系数出现负值，就是稳定性被破坏的直接信号。对 FTCS，中心系数 1 − 2r 在 r > 1/2 时变负，与冯·诺依曼分析结论完全一致。

四、实现要点与常见误区

• 先由 α 和 Δx 算出 Δt 上限，再取一个留有余量的步长（如上限的 0.8 倍），不要贴着临界值运行。
• 边界条件在每一步更新后单独赋值，不要参与内点循环，否则会引入额外误差。
• r = 1/2 时格式恰好稳定但无阻尼余量，实际计算建议 r ≤ 0.4。
• 若必须使用大时间步长，应改用隐式格式（如 Crank–Nicolson），它无条件稳定，但每步需解三对角方程组。