📚 课程目录第2章有限差分法基本思想与网格剖分

2. 有限差分法基本思想与网格剖分

有限差分法求解热方程的核心思路只有一句话：把连续的温度场 T(x,t) 换成有限个离散节点上的温度值，把导数换成差商。一旦完成这一步，偏微分方程就变成了一组代数方程，可以用计算机逐点、逐时刻地推进求解。本节要解决的就是"怎么离散"和"离散后怎么组织计算"这两个前置问题，它们是后续显式格式与隐式格式的共同基础。

一、空间离散：把一根杆切成节点

以一维均匀网格为例。设求解区域为 0 ≤ x ≤ L，取空间步长 Δx，节点编号 i = 0,1,2,…,N，则节点坐标为 x_i = i·Δx，且要求 N·Δx = L。节点 i 处的温度记作 T_i，它代表的是小区间 [x_i−Δx/2, x_i+Δx/2] 上的平均温度近似。i = 0 与 i = N 是边界节点，其温度由边界条件直接给定（第一类边界）或通过边界差分方程求出（第二、三类边界）。

空间步长的选取不是任意的。若杆长 L = 1 m，取 N = 20，则 Δx = 0.05 m；取 N = 100，则 Δx = 0.01 m。Δx 越小，空间分辨率越高，但未知量个数线性增加，计算量也随之上升。工程上常先按"每段内部温度变化是否足够平缓"来估计 Δx，再用网格无关性检验（把 Δx 减半，看结果变化是否小于允许误差）确认。

二、时间离散：把过程切成时间层

时间方向同样离散：取时间步长 Δt，时间层编号 n = 0,1,2,…，时刻 t_n = n·Δt。节点温度写成双下标形式 T_iⁿ，表示第 n 个时间层上第 i 个节点的温度。初始条件给出 T_i⁰（n = 0 整层已知），求解任务就是由已知的第 n 层推出第 n+1 层，一层一层向前推进。

这里有一个关键约定必须提前明确：本文所有差分格式中，时间导数一律用向前差分，空间二阶导数用中心差分，且热扩散系数 α 取正值，热流方向由高温指向低温。参照系与符号约定在后续章节中保持不变。

三、用差商近似导数

由泰勒展开可得三个基本近似式。对时间一阶导数：

∂T/∂t ≈ (T_i^(n+1) − T_i^n) / Δt        （向前差分，截断误差 O(Δt)）

对空间二阶导数，用中心差分精度更高：

∂²T/∂x² ≈ (T_(i−1)^n − 2T_i^n + T_(i+1)^n) / Δx²   （截断误差 O(Δx²)）

把这两式代入一维热方程 ∂T/∂t = α·∂²T/∂x²，并令 r = α·Δt/Δx²（称为网格傅里叶数），就得到显式格式：

T_i^(n+1) = T_i^n + r·(T_(i−1)^n − 2T_i^n + T_(i+1)^n)
          = r·T_(i−1)^n + (1−2r)·T_i^n + r·T_(i+1)^n

这个式子清楚地表明：新时刻某点的温度，只由旧时刻它自己和左右邻居三点决定。若把空间导数取在新时间层 n+1 上，则得到隐式格式，同一时间层上所有未知温度相互耦合，需要解一个三对角线性方程组。两种格式的差别，本质上就是"空间差分取在哪一层"这一个选择。

四、网格与时间步进框架

下图给出了一维均匀网格的节点布局、边界节点位置，以及显式格式中"由第 n 层三点计算第 n+1 层一点"的依赖关系。图中同时标出了 Δx、Δt 与网格傅里叶数 r 的作用位置，便于对照理解离散参数如何进入计算。

五、一个可手算的示例

设一根长 L = 1 m 的均匀杆，热扩散系数 α = 1×10⁻⁴ m²/s，初始温度全为 20 ℃，两端突然分别固定在 100 ℃ 与 0 ℃。取 N = 10，则 Δx = 0.1 m；取 Δt = 50 s，则 r = 1×10⁻⁴ × 50 / 0.1² = 0.5。此时显式格式化为 T_iⁿ⁺¹ = 0.5·T_i−1ⁿ + 0·T_iⁿ + 0.5·T_i+1ⁿ，即新值恰好是左右邻居的平均值。以 i = 1 为例，第一层计算为 0.5×100 + 0.5×20 = 60 ℃；i = 2 为 0.5×20 + 0.5×20 = 20 ℃。逐层推进即可得到温度分布随时间向稳态线性分布演化的过程。这个例子中 r = 0.5 恰好落在显式格式稳定边界 r ≤ 0.5 上，若把 Δt 改为 100 s，r = 1，计算会迅速振荡发散——这正是下一节讨论稳定性时要重点分析的现象。

网格剖分完成后，建议先做两项自检：一是 N·Δx 是否严格等于 L，避免边界节点错位；二是由 α、Δt、Δx 算出的 r 是否落在所选格式的稳定范围内。这两项检查只需几秒钟，却能避免整段计算作废。

六、离散框架的组织方式

把上述内容整理成可执行的框架，包含四个要素：节点数组（长度 N+1，存放当前层温度）、边界处理（每层推进后强制写入边界值）、时间循环（n 从 0 到 n_max，每步由旧层算新层）、结果记录（按需保存若干时间层的温度分布用于后处理）。伪代码如下：

输入：N, L, alpha, dt, t_end, T_left, T_right
dx = L / N;  r = alpha * dt / dx^2
T[0..N] = 初始温度;  T[0] = T_left;  T[N] = T_right
for n = 0 .. t_end/dt − 1:
    Tnew = T                      # 复制旧层
    for i = 1 .. N−1:
        Tnew[i] = r*T[i−1] + (1−2r)*T[i] + r*T[i+1]
    Tnew[0] = T_left;  Tnew[N] = T_right
    T = Tnew                      # 整层更新，不可边算边覆盖
输出：各时刻 T[0..N]

注意伪代码中"先复制、后整层更新"的顺序：如果直接在原数组上逐点覆盖，i+1 点用到的将是本层已更新的新值，格式就悄悄变成了另一种算法。这是初学者最常见的实现错误之一。隐式格式只需把内层循环替换为组装三对角方程组并求解，外层时间推进框架完全不变——这正是本节强调"先搭好网格与时间步进骨架"的意义所在。