3. 二维泊松方程五点差分格式与迭代求解

问题设定与离散格式

在矩形区域 Ω=[0,1]×[0,1] 上求解二维泊松方程 −Δu=f，边界取齐次 Dirichlet 条件 u=0。取均匀网格 x_i=ih、y_j=jh，h=1/N，内点 (i,j) 满足 1≤i,j≤N−1。对拉普拉斯算子用二阶中心差分，得到五点差分格式：

−(u[i+1,j] + u[i−1,j] + u[i,j+1] + u[i,j−1] − 4u[i,j]) / h² = f[i,j]

整理为便于迭代的形式，记 U_i,j 为数值解：

U[i,j] = ( U[i+1,j] + U[i−1,j] + U[i,j+1] + U[i,j−1] + h²·f[i,j] ) / 4

该格式截断误差为 O(h²)，边界值直接取 0，无需额外处理。下面用制造解验证误差阶：取真解 u*(x,y)=sin(πx)sin(πy)，则 f=2π²sin(πx)sin(πy)，边界自然为零。

完整可运行实现

下面用 Python 实现网格生成、右端项计算、两种迭代求解与误差统计。所有参数在本例中定义，可直接运行。

import numpy as np

def solve_poisson(N, tol=1e-10, max_iter=200000):
    h = 1.0 / N
    x = np.linspace(0, 1, N + 1)
    y = np.linspace(0, 1, N + 1)
    X, Y = np.meshgrid(x, y, indexing='ij')
    # 制造解 u* = sin(pi x) sin(pi y)，f = 2 pi^2 sin(pi x) sin(pi y)
    f = 2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
    u_exact = np.sin(np.pi * X) * np.sin(np.pi * Y)

    # 初始化内点为 0，边界保持 0
    U = np.zeros((N + 1, N + 1))
    h2f = h * h * f

    # Jacobi
    Uj = U.copy()
    it_j = 0
    for it_j in range(1, max_iter + 1):
        Un = Uj.copy()
        Un[1:-1, 1:-1] = 0.25 * (Uj[2:, 1:-1] + Uj[:-2, 1:-1]
                                 + Uj[1:-1, 2:] + Uj[1:-1, :-2]
                                 + h2f[1:-1, 1:-1])
        diff = np.max(np.abs(Un - Uj))
        Uj = Un
        if diff < tol:
            break

    # Gauss-Seidel（就地更新）
    Ug = U.copy()
    it_g = 0
    for it_g in range(1, max_iter + 1):
        diff = 0.0
        for i in range(1, N):
            for j in range(1, N):
                old = Ug[i, j]
                Ug[i, j] = 0.25 * (Ug[i+1, j] + Ug[i-1, j]
                                   + Ug[i, j+1] + Ug[i, j-1]
                                   + h2f[i, j])
                diff = max(diff, abs(Ug[i, j] - old))
        if diff < tol:
            break

    err_j = np.max(np.abs(Uj - u_exact))
    err_g = np.max(np.abs(Ug - u_exact))
    return it_j, it_g, err_j, err_g

for N in [8, 16, 32, 64]:
    it_j, it_g, ej, eg = solve_poisson(N)
    print(f"N={N:3d}  Jacobi迭代={it_j:6d}  GS迭代={it_g:6d}  "
          f"Jacobi误差={ej:.3e}  GS误差={eg:.3e}")

结果判读与误差阶验证

运行后应观察到两类现象。第一，迭代次数：Gauss-Seidel 的迭代次数约为 Jacobi 的一半，与谱半径 ρ_GS=ρ_J² 的结论一致。第二，误差阶：由于两种迭代都收敛到同一离散解，最大误差应随 N 加倍而约缩小为 1/4，即误差比 ≈4，对应 O(h²)。