4. 二维泊松方程九点差分格式精度提升实践
问题设定与制造解
在单位正方形 Ω=[0,1]² 上求解二维泊松方程 −Δu=f，取制造解 u(x,y)=sin(πx)sin(πy)，则右端项 f=2π²sin(πx)sin(πy)，边界上 u=0。该解在边界附近曲率变化剧烈，正好用来检验高阶格式对边界附近精度的改善。
为便于比较，两种格式使用相同网格：N×N 内部节点，h=1/(N+1)，节点坐标 xi=ih，yj=jh。误差用离散最大范数 E∞=max|uij−u(xi,yj)| 与离散 L² 范数衡量，收敛阶用相邻网格误差比的对数估计：p≈log(Eh/Eh/2)/log2。
五点格式与九点紧致格式的构造
五点格式直接对 −(uxx+uyy) 作中心二阶差分：
(-Δ_h^5 u)_ij = (4u_ij - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1}) / h²其截断误差为 O(h²)。九点紧致格式利用 uxx、uyy 的四阶紧致算子组合，并借助方程本身消去高阶导数项。记 δx²uij=ui+1,j−2uij+ui−1,j，δy² 同理，则九点格式可写为：
(1 + h²/12 · (δ_x² + δ_y²)/h²) 作用在 (-Δ_h^5 u) 上 = f_ij + h²/12 · (δ_x² + δ_y²)f_ij / h²展开后每个内部节点涉及 3×3 共九个点（中心、四邻、四对角），故称九点格式。其截断误差为 O(h⁴)，但边界附近需要特殊处理：由于格式模板在距边界一层的节点上会引用域外点，实践中对最靠近边界的节点仍用五点格式，内部节点用九点格式，这种"边界降阶"是紧致格式的典型工程处理方式。
可运行的实现代码
下面给出完整的 Python 实现，包含制造解、右端项、两种格式的系数矩阵组装、直接求解（小网格用稠密求解便于验证）以及误差与收敛阶计算。代码自洽，可直接运行。
import numpy as np
def exact(x, y):
return np.sin(np.pi * x) * np.sin(np.pi * y)
def rhs(x, y):
return 2.0 * np.pi**2 * np.sin(np.pi * x) * np.sin(np.pi * y)
def build_grid(N):
h = 1.0 / (N + 1)
x = np.linspace(h, 1 - h, N)
y = np.linspace(h, 1 - h, N)
X, Y = np.meshgrid(x, y, indexing='ij')
return h, X, Y
def solve_five_point(N):
h, X, Y = build_grid(N)
n = N * N
A = np.zeros((n, n))
b = np.zeros(n)
def idx(i, j):
return i * N + j
for i in range(N):
for j in range(N):
k = idx(i, j)
A[k, k] = 4.0 / h**2
b[k] = rhs(X[i, j], Y[i, j])
if i > 0: A[k, idx(i-1, j)] = -1.0 / h**2
if i < N-1: A[k, idx(i+1, j)] = -1.0 / h**2
if j > 0: A[k, idx(i, j-1)] = -1.0 / h**2
if j < N-1: A[k, idx(i, j+1)] = -1.0 / h**2
u = np.linalg.solve(A, b).reshape(N, N)
return h, X, Y, u
def solve_nine_point(N):
h, X, Y = build_grid(N)
n = N * N
A = np.zeros((n, n))
b = np.zeros(n)
def idx(i, j):
return i * N + j
# 九点格式系数（乘以 h² 归一化后）：
# 中心 -20/6，四邻 4/6，四对角 1/6，右端为 h²*(f + h²/12 * Δf)
# 这里采用等价形式：对 -Δu=f 两边作用 (1 + h²/12 Δ) 得
# (-Δ_h^5 u) + h²/12 * (-Δ_h^5)(δ_x²+δ_y²)/h² u = f + h²/12 Δf
# 展开后中心 -20/6，四邻 4/6，四对角 1/6（均除以 h²）
c0, c1, c2 = -20.0/6.0, 4.0/6.0, 1.0/6.0
for i in range(N):
for j in range(N):
k = idx(i, j)
# 边界层（距边界一层）用五点格式
boundary_layer = (i == 0 or i == N-1 or j == 0 or j == N-1)
if boundary_layer:
A[k, k] = 4.0 / h**2
b[k] = rhs(X[i, j], Y[i, j])
if i > 0: A[k, idx(i-1, j)] = -1.0 / h**2
if i < N-1: A[k, idx(i+1, j)] = -1.0 / h**2
if j > 0: A[k, idx(i, j-1)] = -1.0 / h**2
if j < N-1: A[k, idx(i, j+1)] = -1.0 / h**2
else:
A[k, k] = c0 / h**2
for di, dj, c in [(-1,0,c1),(1,0,c1),(0,-1,c1),(0,1,c1),
(-1,-1,c2),(-1,1,c2),(1,-1,c2),(1,1,c2)]:
A[k, idx(i+di, j+dj)] = c / h**2
# 右端：f + h²/12 * Δf，Δf 用中心差分
f0 = rhs(X[i, j], Y[i, j])
fxx = (rhs(X[i+1,j],Y[i+1,j]) - 2*f0 + rhs(X[i-1,j],Y[i-1,j])) / h**2
fyy = (rhs(X[i,j+1],Y[i,j+1]) - 2*f0 + rhs(X[i,j-1],Y[i,j-1])) / h**2
b[k] = f0 + h**2/12.0 * (fxx + fyy)
u = np.linalg.solve(A, b).reshape(N, N)
return h, X, Y, u
def error_norm(u, X, Y):
ue = exact(X, Y)
return np.max(np.abs(u - ue)), np.sqrt(np.mean((u - ue)**2))
if __name__ == "__main__":
print(f"{'N':>4} {'五点E_inf':>12} {'五点p':>8} {'九点E_inf':>12} {'九点p':>8}")
prev5 = prev9 = None
for N in [8, 16, 32, 64]:
h5, X5, Y5, u5 = solve_five_point(N)
h9, X9, Y9, u9 = solve_nine_point(N)
e5, _ = error_norm(u5, X5, Y5)
e9, _ = error_norm(u9, X9, Y9)
p5 = np.log2(prev5/e5) if prev5 else float('nan')
p9 = np.log2(prev9/e9) if prev9 else float('nan')
print(f"{N:>4} {e5:>12.3e} {p5:>8.2f} {e9:>12.3e} {p9:>8.2f}")
prev5, prev9 = e5, e9
结果判读与精度分析
运行上述代码（示例网格 N=8,16,32,64）可观察到以下典型现象，这些是由格式截断误差直接推出的预期结果，而非实测数据：
- 五点格式：E∞ 随 h 减小按约 2 阶下降，相邻网格误差比接近 4，估计收敛阶 p≈2。
- 九点紧致格式：内部节点误差按约 4 阶下降，但边界层降阶为五点后，整体收敛阶通常介于 2 与 4 之间；在 N 较大时，边界层误差占比下降，整体阶数趋近 3～4。这正是"边界附近精度"成为紧致格式工程应用关键的原因。
- 相同网格对比：在 N=32 时，九点格式的 E∞ 通常比五点格式小 1～2 个数量级（取决于边界层处理方式），且误差分布更集中于边界附近。
边界附近精度的可视化判读
将误差按节点位置绘制热力图可直观看到：五点格式的误差在边界附近与内部量级相近，而九点格式的误差在内部极小、在边界一层明显抬升。若把边界层也升级为高阶处理，误差分布会趋于均匀。判断格式是否"真正高阶"的关键，不是看内部最大误差，而是看误差随 h 的整体衰减斜率是否与理论阶一致。