4. 有限差分法离散一维泊松方程

一维泊松方程的标准形式为

−u″(x) = f(x)，  x ∈ (0, L)，  u(0) = a，  u(L) = b

其中 f(x) 是已知源项，u(x) 是待求函数，a、b 为两端狄利克雷边界值。本节围绕“如何把连续的微分方程变成可计算的代数方程组”展开：先推导三点差分格式，再分析截断误差，然后处理边界条件，最后用一个可手算的算例走通全流程。全文约定 x 轴向右为正，网格节点从左到右编号 i = 0, 1, …, N，步长 h = L/N，节点坐标 x_i = i·h。

一、三点差分格式的推导

在内部节点 x_i 处，把 u(x_i±h) 分别做泰勒展开：

u(x_i+h) = u_i + h·u′_i + h²/2·u″_i + h³/6·u‴_i + h⁴/24·u⁗_i + …

u(x_i−h) = u_i − h·u′_i + h²/2·u″_i − h³/6·u‴_i + h⁴/24·u⁗_i − …

两式相加，一阶导数项与三阶导数项相消，得到二阶中心差分：

u″_i = (u_i−1 − 2u_i + u_i+1) / h² + O(h²)

代入 −u″ = f，得到内部节点的三点差分方程：

−u_i−1 + 2u_i − u_i+1 = h²·f_i，  i = 1, 2, …, N−1

每个内部节点只与左右相邻节点耦合，因此 N−1 个未知量 u₁…u_N−1 构成一个三对角线性方程组。写成矩阵形式为 A·U = F，其中 A 的主对角线全为 2，次对角线与超对角线全为 −1，右端项 F_i = h²·f_i，边界值 u₀ = a、u_N = b 移到右端后并入 F₁ 与 F_N−1。该矩阵严格对角占优，方程组有唯一解，可用追赶法（Thomas 算法）以 O(N) 复杂度求解。

二、截断误差分析

把上面两式相减可得一阶导数中心差分，而二阶导数的截断误差需保留到四阶项。将两式相加后除以 h²：

(u_i−1 − 2u_i + u_i+1)/h² = u″_i + h²/12·u⁗_i + O(h⁴)

因此局部截断误差 τ_i = h²/12·u⁗(ξ_i)，格式具有二阶精度。这带来两个实用结论：其一，步长减半时误差约缩小为原来的四分之一，可用“网格加密误差比接近 4”来验证实现是否正确；其二，误差与解的四阶导数成正比，若源项 f 变化剧烈或解曲率大，粗网格上的误差会明显偏大，需要更细的网格或更高阶格式。

三、边界条件的处理

狄利克雷边界最简单：u₀、u_N 已知，直接移入右端项，不参与求解。例如 i = 1 的方程变为 2u₁ − u₂ = h²f₁ + a，i = N−1 的方程变为 −u_N−2 + 2u_N−1 = h²f_N−1 + b。

若边界给定的是导数（诺伊曼条件），则需引入虚拟节点。以左端 u′(0) = g 为例，在 x₀ 处用中心差分 u′(0) ≈ (u₁ − u₋₁)/(2h) = g，得 u₋₁ = u₁ − 2hg。再对 i = 0 写三点方程 −u₋₁ + 2u₀ − u₁ = h²f₀，消去 u₋₁ 后得到 2u₀ − 2u₁ = h²f₀ − 2hg。这样边界节点也进入方程组，矩阵仍保持三对角结构，只是首行元素变为 (2, −2)。这种虚拟节点法保持二阶精度，优于把边界导数用单侧一阶差分硬代入的做法。

四、简单算例

取 L = 1，f(x) = 1，边界 u(0) = 0、u(1) = 0。该问题的解析解为 u(x) = x(1−x)/2，便于对照误差。取 N = 4，则 h = 0.25，h² = 0.0625，内部节点 i = 1, 2, 3，右端项均为 0.0625。

方程组为：

2u1 − u2        = 0.0625 + 0 = 0.0625
−u1 + 2u2 − u3  = 0.0625
       −u2 + 2u3 = 0.0625 + 0 = 0.0625

由对称性 u₁ = u₃，代入解得 u₂ = 0.09375，u₁ = u₃ = 0.078125。与解析解对比：x = 0.25 处精确值 0.09375，x = 0.5 处精确值 0.125。可见节点值存在偏差，最大误差出现在中点，约 0.03125。若把 N 加密到 8，h² 变为 0.015625，误差应约降为原来的四分之一，这正是二阶收敛的体现。

从算例可以看出，离散求解的完整链条是：确定网格与步长 → 组装三对角矩阵与右端项 → 处理边界 → 用追赶法求解 → 与解析解或加密网格结果对比验证。实际编程时，右端项务必乘上 h²，边界值务必正确并入首末方程，这两处是最常见的错误来源。截断误差分析则告诉我们：当结果与参考解偏差较大时，优先检查网格是否足够细，而不是盲目怀疑算法本身。