4. 一维泊松方程完整求解流程

本节把前面分散的知识点串成一条可执行的完整链路：给定一维泊松方程，从定义网格开始，到构造函数空间、施加边界条件、写出变分形式、调用 solve，最后把数值解与解析解放在一起比较，用误差指标判断精度是否达标。为了让流程可复现，我们采用一个解析解已知的算例，这样每一步的结果都能被独立检验。

算例设定与解析解

考虑区间 Ω = [0, 1] 上的一维泊松方程：

-u''(x) = f(x),  x ∈ (0, 1)
u(0) = 0,  u(1) = 0

取真解 u(x) = sin(πx)，则源项 f(x) = π² sin(πx)。这个算例的好处是：真解在端点为零，与齐次 Dirichlet 边界条件天然匹配；同时 sin(πx) 光滑，有限元解应当以最优阶收敛，便于我们观察误差随网格加密的下降趋势。注意符号约定：本文统一采用 -u'' = f 的形式，弱形式中会出现正号的双线性项，若你习惯 u'' = f 的写法，所有符号需整体取反，切勿混用。

第一步：定义网格

一维问题用 UnitIntervalMesh 生成区间网格，参数是区间被划分的单元数。先取较粗的 8 个单元做流程验证，再加密到 32、128 观察收敛：

from dolfin import *
import numpy as np

N = 8
mesh = UnitIntervalMesh(N)
print("顶点数:", mesh.num_vertices(), " 单元数:", mesh.num_cells())

网格对象本身只描述几何剖分，不携带任何未知量信息。一维情况下每个单元是线段，顶点即自由度载体，这一点在后续判断自由度数量时很有用。

第二步：构造函数空间

泊松方程是二阶方程，弱形式中只出现一阶导数，因此只需 C⁰ 连续的分片线性函数即可，选用一阶 Lagrange 单元：

V = FunctionSpace(mesh, "CG", 1)

对 N = 8 的网格，区间被 8 个单元分成 9 个顶点，故 V.dim() 应为 9。这个数字值得在代码里打印核对一次，它是后面组装矩阵维度和施加边界条件时的重要参照。

第三步：边界条件

两个端点都是本质边界条件，用 DirichletBC 分别指定。一维中 on_boundary 会同时命中 x=0 和 x=1，因此必须用坐标判断把两端区分开：

u_D = Constant(0.0)
bc_left  = DirichletBC(V, u_D, "on_boundary && near(x[0], 0.0)")
bc_right = DirichletBC(V, u_D, "on_boundary && near(x[0], 1.0)")
bcs = [bc_left, bc_right]

若只写 on_boundary 而不加坐标条件，两端都会被赋同一个值——本例恰好都是 0 所以看不出问题，但换成非齐次边界就会静默出错。这是初学者最常见的陷阱之一，建议养成显式区分边界的习惯。

第四步：变分形式

用测试函数 v 乘方程两端并在 Ω 上积分，对 -u''v 分部积分，利用 v 在边界为零消去边界项，得到弱形式：求 u ∈ V 且满足边界条件，使得对所有 v ∈ V 有

∫ u' v' dx = ∫ f v dx

对应到代码，dx 是一维体积测度，grad 在一维下退化为对 x 的导数：

u = TrialFunction(V)
v = TestFunction(V)
f = Expression("pi*pi*sin(pi*x[0])", pi=np.pi, degree=4)

a = inner(grad(u), grad(v)) * dx
L = f * v * dx

源项用 Expression 给出，degree=4 指定插值精度，保证积分时对 f 的数值求积足够准确，避免把求积误差误判为离散误差。

第五步：调用 solve 并可视化

把双线性形式、线性形式、未知函数和边界条件交给 solve，它会自动完成组装、施加边界条件与线性求解：

uh = Function(V)
solve(a == L, uh, bcs)

# 与解析解对比
u_exact = Expression("sin(pi*x[0])", pi=np.pi, degree=4)
err = errornorm(u_exact, uh, "L2")
print("L2 误差:", err)

可视化方面，一维解最直观的方式是采样后画曲线。把网格顶点坐标与对应自由度值取出，与解析解在同一组点上比较：

xs = np.linspace(0, 1, 200)
uh_vals = np.array([uh(Point(x)) for x in xs])
ex_vals = np.sin(np.pi * xs)

若在 Jupyter 中，也可直接用 plot(uh) 快速查看形状；但要做定量比较，采样成数组更可靠。

上图刻意画的是 N=8 的粗网格结果：蓝色折线是分片线性有限元解，红色虚线是解析解。可以看到折线在峰值附近明显偏低，这正是离散误差的直观体现。随着 N 增大，折线会逐步贴合曲线，误差以二阶速率下降。

精度验证：收敛阶的定量判断

单次误差数值本身意义有限，真正有说服力的是收敛阶。对一阶 Lagrange 单元，L2 误差理论上满足 ‖u - u_h‖₂ ≤ C h²，其中 h = 1/N。把 N 依次取 8、16、32、64，计算相邻两次误差的比值，若比值接近 4，说明收敛阶约为 2，实现正确：

for N in [8, 16, 32, 64]:
    mesh = UnitIntervalMesh(N)
    V = FunctionSpace(mesh, "CG", 1)
    ...  # 重复上述组装与求解
    e = errornorm(u_exact, uh, "L2")
    print(N, e)

一个合理的示例结果序列大致是：N=8 时误差约 1e-2 量级，N=16 约 2.5e-3，N=32 约 6e-4，N=64 约 1.5e-4，每次缩小约 4 倍。若你实际运行得到的比值明显偏离 4——比如只缩小 2 倍——通常意味着源项插值精度不足、边界条件漏加，或变分形式符号写反，应逐项回查。

至此，一条完整的一维泊松方程求解链路就闭合了：网格决定离散精度，函数空间决定逼近能力，边界条件决定问题适定性，变分形式决定方程等价性，solve 完成数值实现，而误差与收敛阶则给出“结果是否可信”的最终判据。把这套流程固化成脚本，换一个源项或边界条件只需改动少数几行，这正是有限元框架的价值所在。