用稳态求解器计算非齐次边界的 Helmholtz 方程

本文根据 SciML 的 FiniteVolumeMethod.jl 官方教程整理并中文化。教程页未署个人作者；候选记录将作者记为“FiniteVolumeMethod.jl 文档维护者（页面未署个人名）”。

• 来源：https://docs.sciml.ai/FiniteVolumeMethod/stable/tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions
• 原文标题：Helmholtz Equation with Inhomogeneous Boundary Conditions
• 分类：数据处理
• 版本与环境：来源地址为 stable 文档。页面注明使用 Documenter.jl 1.8.0、Julia 1.11.2 于 2025 年 1 月 1 日生成；这不是 FiniteVolumeMethod.jl 的版本号，也不代表文章首发日期。教程没有锁定 FiniteVolumeMethod 或其他依赖的具体包版本。
• 许可说明：页面注明由 Literate.jl 生成，正文未列出文章转载许可。
• 原文图示：教程将计算结果绘制为三角网格上的等值填色图；此处保留代码与文字说明，未附图。

问题与边界条件

在方形区域 [-1, 1]² 上求 Helmholtz 方程的稳态解：

∇²u(x) + u(x) = 0，x ∈ [-1, 1]²

边界条件为沿外法向的导数等于 1：

∂u/∂n = 1，x ∈ ∂([-1, 1]²)

这不是齐次 Neumann 边界，因为边界导数不为零。用有限体积方法建模时，首先需要把它换算成该库使用的通量表达。

创建三角网格

教程使用 DelaunayTriangulation.jl 与 FiniteVolumeMethod.jl，在正方形区域上创建 125 × 125 的矩形剖分，并通过 FVMGeometry 构造有限体积几何：

using DelaunayTriangulation, FiniteVolumeMethod

tri = triangulate_rectangle(-1, 1, -1, 1, 125, 125, single_boundary=true)
mesh = FVMGeometry(tri)

页面输出网格包含 15,625 个控制体积、30,752 个三角形和 46,376 条边。随后定义 Neumann 边界条件：

BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> -one(u), Neumann)

处理通量的符号

原方程给定的是法向导数：

∇u · n = 1

而本例的有限体积通量定义为 q = -∇u。因此，将原条件改写为通量形式时：

q · n = -1

所以回调函数返回的是 -one(u)，而不是 +one(u)。这里的负号来自通量定义；若直接把导数条件数值 1 传给 Neumann 边界，会弄反边界通量。

把方程包装成稳态问题

FiniteVolumeMethod 的时变形式可写为：

∂u/∂t + ∇·q = S

或等价地：

∂u/∂t = ∇·[D∇u] + S

稳态形式则为：

∇·q = S

或：

∇·[D∇u] + S = 0

对于本例，扩散系数 D = 1，源项 S = u。先定义扩散项、源项、全零初始估计，再设置无穷终止时间并创建 FVMProblem。稳态问题中的 initial_condition 是供非线性求解器使用的初始猜测；final_time 则设为 Inf。最后用 SteadyFVMProblem 包装该有限体积问题：

diffusion_function = (x, y, t, u, p) -> one(u)
source_function = (x, y, t, u, p) -> u
initial_condition = zeros(DelaunayTriangulation.num_solid_vertices(tri))
final_time = Inf

prob = FVMProblem(mesh, BCs;
    diffusion_function,
    source_function,
    initial_condition,
    final_time)

steady_prob = SteadyFVMProblem(prob)

求解：先 Newton，再决定是否细化

教程先用 NonlinearSolve.jl 的 NewtonRaphson() 求解稳态问题：

using NonlinearSolve

sol = solve(steady_prob, NewtonRaphson())

如果希望再用稳态差分方法细化，可以将 Newton 得到的解复制到原问题的初始猜测中，然后调用 SteadyStateDiffEq.jl 的 DynamicSS，并使用 TRBDF2 时间步进器与 KLU 稀疏线性求解器：

copyto!(prob.initial_condition, sol.u)
using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq

sol = solve(
    steady_prob,
    DynamicSS(TRBDF2(linsolve=KLUFactorization()))
)

教程也提及一种可选写法，可给 DynamicSS 指定相对容差：

DynamicSS(TRBDF2(linsolve=KLUFactorization()), reltol=1e-4)

本例中 Newton 的结果被复制回 prob.initial_condition；由于稳态包装对象使用同一个初始条件，更新后可作为后续求解的起点。教程展示 DynamicSS 的二次求解，但也明确说明，对于这个问题似乎不需要这一步修正。因此，应把它看作可选的进一步求稳态方式，而不是每个问题都必须套用的固定流程。

页面展示的求解结果报告 retcode: Success，解向量有 15,625 个 Float64 元素。原文仅显示向量首尾若干项，例如开头为 -1.284086882618605、-1.3001602399337975、-1.3160471807169647、-1.3317536108619856、-1.347278372873141；末尾对称地回到约 -1.284086882618604。这里转述的是文档页面所列输出，不是本次运行或独立精度验证的结果。

绘制解的等值图

教程使用 CairoMakie.jl 的 tricontourf 将三角网格上的解画成填色等值图，色阶范围为 -2.5 到 -1.0，步长 0.15，色图为 matter：

using CairoMakie

fig, ax, sc = tricontourf(
    tri,
    sol.u,
    levels=-2.5:0.15:-1.0,
    colormap=:matter
)
fig

输出图展示离散解在整个方形区域内的空间分布。本文保留生成图的代码，不复制原文图片。

完整代码

下面合并教程中的核心步骤，保留其求解顺序与参数：

using DelaunayTriangulation, FiniteVolumeMethod

tri = triangulate_rectangle(-1, 1, -1, 1, 125, 125, single_boundary=true)
mesh = FVMGeometry(tri)

BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> -one(u), Neumann)

diffusion_function = (x, y, t, u, p) -> one(u)
source_function = (x, y, t, u, p) -> u
initial_condition = zeros(DelaunayTriangulation.num_solid_vertices(tri))
final_time = Inf

prob = FVMProblem(mesh, BCs;
    diffusion_function,
    source_function,
    initial_condition,
    final_time)

steady_prob = SteadyFVMProblem(prob)

using NonlinearSolve
sol = solve(steady_prob, NewtonRaphson())
copyto!(prob.initial_condition, sol.u) # this also changes steady_prob's initial condition
using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve=KLUFactorization())))

using CairoMakie
fig, ax, sc = tricontourf(tri, sol.u, levels=-2.5:0.15:-1.0, colormap=:matter)
fig

实践要点

• 转换 Neumann 条件时，要先确认库的通量符号定义。本例 q = -∇u，因此法向导数为 1 对应边界通量 -1。
• SteadyFVMProblem 包装的是稳态计算；initial_condition 此时提供非线性求解初值，final_time 设为 Inf。
• 本例把扩散系数写成 1、源项写成 u。方程要与该库使用的稳态形式对应起来。
• 可以先用 NewtonRaphson() 获得解，再把它传给 DynamicSS 作为进一步稳态求解的起点；教程指出，本例不一定需要第二阶段。
• 等值图由三角网格、解向量、色阶范围和色图共同生成。

来源信息

• 作者：FiniteVolumeMethod.jl 文档维护者（原文页面未标个人姓名）
• 原文标题：Helmholtz Equation with Inhomogeneous Boundary Conditions
• 文档构建信息：Documenter.jl 1.8.0，Julia 1.11.2，构建时间为 2025-01-01
• 原文地址：https://docs.sciml.ai/FiniteVolumeMethod/stable/tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions
• 文章未注明 FiniteVolumeMethod.jl 的包版本；构建时间不作为发表日期