1. 一维稳态扩散方程有限差分离散与验证

问题设定与离散格式

考虑一维稳态扩散方程（取扩散系数 \(D=1\)）：

\(-\dfrac{d^2u}{dx^2}=f(x),\quad x\in(0,1)\)

边界条件取 Dirichlet 型：\(u(0)=0,\ u(1)=0\)。在均匀网格 \(x_i=ih,\ i=0,\dots,N,\ h=1/N\) 上，对内部节点 \(i=1,\dots,N-1\) 用二阶中心差分逼近二阶导数：

\(-\dfrac{u_{i-1}-2u_i+u_{i+1}}{h^2}=f_i\)

整理得线性方程组：

\(-u_{i-1}+2u_i-u_{i+1}=h^2 f_i\)

边界节点 \(u_0=0,\ u_N=0\) 直接代入，不参与未知量求解。这样得到 \(N-1\) 阶三对角系统，可用追赶法（Thomas 算法）\(O(N)\) 求解。

制造解验证思路

为验证离散格式的收敛阶，选取已知解析解 \(u^*(x)=\sin(\pi x)\)，反推源项：

\(f(x)=-u^{*\prime\prime}(x)=\pi^2\sin(\pi x)\)

数值解 \(u_h\) 与解析解在网格点上的最大误差记为 \(E(h)=\max_i|u_i-u^*(x_i)|\)。对二阶格式，理论上 \(E(h)=Ch^2\)，即网格加密一倍误差约降为 1/4。用两组网格 \(N\) 与 \(2N\) 计算收敛阶：

\(p\approx\dfrac{\log(E(h)/E(h/2))}{\log 2}\)

关键实现（Python）

下面给出可直接运行的完整实现：构造三对角系统、用 Thomas 算法求解，并在两组网格上计算误差与收敛阶。

import numpy as np

def solve_diffusion(N):
    """求解 -u'' = pi^2 sin(pi x), u(0)=u(1)=0，返回网格点与数值解"""
    h = 1.0 / N
    x = np.linspace(0.0, 1.0, N + 1)
    # 内部未知量 u_1..u_{N-1}
    n = N - 1
    f = np.pi**2 * np.sin(np.pi * x[1:N])
    b = h**2 * f

    # 三对角：下对角 a、主对角 d、上对角 c
    a = -np.ones(n - 1)
    d =  2.0 * np.ones(n)
    c = -np.ones(n - 1)

    # Thomas 算法
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / d[0] if n > 1 else 0.0
    dp[0] = b[0] / d[0]
    for i in range(1, n):
        m = d[i] - a[i - 1] * cp[i - 1]
        dp[i] = (b[i] - a[i - 1] * dp[i - 1]) / m
        if i < n - 1:
            cp[i] = c[i] / m
    u_int = np.zeros(n)
    u_int[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        u_int[i] = dp[i] - cp[i] * u_int[i + 1]

    u = np.zeros(N + 1)
    u[1:N] = u_int
    return x, u

def max_error(N):
    x, u = solve_diffusion(N)
    u_exact = np.sin(np.pi * x)
    return np.max(np.abs(u - u_exact))

for N in [10, 20, 40, 80]:
    e = max_error(N)
    print(f"N={N:3d}  h={1/N:.4f}  max_err={e:.3e}")

# 收敛阶
e1 = max_error(20); e2 = max_error(40)
p = np.log(e1 / e2) / np.log(2)
print(f"估计收敛阶 p ≈ {p:.3f}")

结果判读与边界条件影响

预期输出（示例数据，由上述代码直接计算得到）：

误差比接近 4、收敛阶 \(p\approx 2\)，说明中心差分格式在光滑解上达到二阶精度。