用 FVMSystem 联立求解 Brusselator 反应扩散方程-未完纪

原文标题： A Reaction-Diffusion Brusselator System of PDEs 作者： FiniteVolumeMethod.jl 文档维护者（原页面未署个人名） 来源： SciML / FiniteVolumeMethod.jl 官方文档 原文： https://docs.sciml.ai/FiniteVolumeMethod/stable/tutorials/reaction_diffusion_brusselator_system_of_pdes 文档条件： Notion 核验记录标注读取的是 stable 文档；页面设置使用 Julia 1.11.2，Documenter.jl 1.8.0 于 2025-01-01 生成。页面没有锁定 FiniteVolumeMethod 或其他依赖的具体版本，也没有注明文章转载许可。

本教程展示如何用 FVMSystem 求解两个相互耦合的偏微分方程（PDE），并把不同变量的通量边界、反应源项和初始条件分别定义后组合起来。示例在正方形区域上求解 Brusselator 反应扩散系统。

方程和解析解

定义区域为 \((x,y)\in[0,1]^2\)，方程为：

\[ \begin{aligned} \frac{\partial \Phi}{\partial t} &= \frac14\nabla^2\Phi+\Phi^2\Psi-2\Phi,\\ \frac{\partial \Psi}{\partial t} &= \frac14\nabla^2\Psi-\Phi^2\Psi+\Phi. \end{aligned} \]

为了直接构造初值和边界条件，示例选用已知解析形式的解：

\[ \Phi(x,y,t)=e^{-x-y-t/2},\qquad \Psi(x,y,t)=e^{x+y+t/2}. \]

因此初始条件是：

\[ \Phi(x,y,0)=e^{-x-y},\qquad \Psi(x,y,0)=e^{x+y}. \]

教程在正方形四条边上混用 Dirichlet 与 Neumann 条件：

边界	\(\Phi\) 条件	\(\Psi\) 条件
\(y=0\)	\(\partial_y\Phi=-e^{-x-t/2}\)	\(\Psi=e^{x+t/2}\)
\(x=1\)	\(\partial_x\Phi=-e^{-1-y-t/2}\)	\(\partial_x\Psi=e^{1+y+t/2}\)
\(y=1\)	\(\Phi=e^{-1-x-t/2}\)	\(\partial_y\Psi=e^{1+x+t/2}\)
\(x=0\)	\(\partial_x\Phi=-e^{-y-t/2}\)	\(\Psi=e^{y+t/2}\)

有限体积法通过法向通量表达 Neumann 条件。对 \(\Phi\)、\(\Psi\) 分别记通量为 \(\mathbf q_1\)、\(\mathbf q_2\)，边界需要写成 \(\mathbf q_i\cdot\mathbf n=f(x,y,t)\) 的形式，其中 \(\mathbf n\) 是外法向。因而不能只照搬单变量教程中的反应扩散简写：本文档版本对多变量系统要求使用守恒通量形式，这是其 API 的限制之一。下面的回调按底、右、顶、左的顺序填写，并已把外法向带来的符号包含在通量值中。

建立网格

先载入有限体积方法和三角剖分包，在单位正方形上创建网格。原示例使用 100×100 的矩形划分：


using FiniteVolumeMethod, DelaunayTriangulation

tri = triangulate_rectangle(0, 1, 0, 1, 100, 100, single_boundary=false)
mesh = FVMGeometry(tri)

页面展示此网格包含 10,000 个控制体积、19,602 个三角形和 29,601 条边。以上为原文示例输出，不表示本任务实际运行过代码。

分别定义两个变量的边界条件

对 PDE 系统，每个变量都要建立自己的边界条件。回调仍接收 (x, y, t, u, p)，但 u 是各变量解值组成的向量或元组，而不是单个数值。本例的边界函数不需要读取 u，但定义其他耦合边界时要留意此签名变化。


Φ_bot = (x, y, t, u, p) -> -1 / 4 * exp(-x - t / 2)
Φ_right = (x, y, t, u, p) -> 1 / 4 * exp(-1 - y - t / 2)
Φ_top = (x, y, t, u, p) -> exp(-1 - x - t / 2)
Φ_left = (x, y, t, u, p) -> -1 / 4 * exp(-y - t / 2)
Φ_bc_fncs = (Φ_bot, Φ_right, Φ_top, Φ_left)
Φ_bc_types = (Neumann, Neumann, Dirichlet, Neumann)
Φ_BCs = BoundaryConditions(mesh, Φ_bc_fncs, Φ_bc_types)

对应的类型顺序为底边 Neumann、右边 Neumann、顶边 Dirichlet、左边 Neumann。Ψ 的边界数据如下：


Ψ_bot = (x, y, t, u, p) -> exp(x + t / 2)
Ψ_right = (x, y, t, u, p) -> -1 / 4 * exp(1 + y + t / 2)
Ψ_top = (x, y, t, u, p) -> -1 / 4 * exp(1 + x + t / 2)
Ψ_left = (x, y, t, u, p) -> exp(y + t / 2)
Ψ_bc_fncs = (Ψ_bot, Ψ_right, Ψ_top, Ψ_left)
Ψ_bc_types = (Dirichlet, Neumann, Neumann, Dirichlet)
Ψ_BCs = BoundaryConditions(mesh, Ψ_bc_fncs, Ψ_bc_types)

即底边 Dirichlet、右边 Neumann、顶边 Neumann、左边 Dirichlet。Neumann 回调返回的数值是外法向通量，而不是未经变换的坐标方向导数。

定义通量、反应项和初始值

通量回调中的梯度系数 α、β、γ 在系统问题里改为按变量排列的元组。索引 1 对应 Φ，索引 2 对应 Ψ。源项收到的状态也以 (Φ, Ψ) 元组给出：


Φ_q = (x, y, t, α, β, γ, p) -> (-α[1] / 4, -β[1] / 4)
Ψ_q = (x, y, t, α, β, γ, p) -> (-α[2] / 4, -β[2] / 4)
Φ_S = (x, y, t, (Φ, Ψ), p) -> Φ^2 * Ψ - 2Φ
Ψ_S = (x, y, t, (Φ, Ψ), p) -> -Φ^2 * Ψ + Φ

使用解析解生成网格节点上的初始值：


Φ_exact = (x, y, t) -> exp(-x - y - t / 2)
Ψ_exact = (x, y, t) -> exp(x + y + t / 2)
Φ₀ = [Φ_exact(x, y, 0) for (x, y) in DelaunayTriangulation.each_point(tri)]
Ψ₀ = [Ψ_exact(x, y, 0) for (x, y) in DelaunayTriangulation.each_point(tri)];

创建问题并联立求解

先为每个变量分别构造 FVMProblem，再将二者传入 FVMSystem：


Φ_prob = FVMProblem(mesh, Φ_BCs; flux_function=Φ_q, source_function=Φ_S,
    initial_condition=Φ₀, final_time=5.0)

Ψ_prob = FVMProblem(mesh, Ψ_BCs; flux_function=Ψ_q, source_function=Ψ_S,
    initial_condition=Ψ₀, final_time=5.0)

system = FVMSystem(Φ_prob, Ψ_prob)

随后通过 OrdinaryDiffEq 和 LinearSolve 求解，使用 TRBDF2 时间积分算法、KLU 稀疏矩阵分解，并每隔 1 个时间单位保存一次结果：


using OrdinaryDiffEq, LinearSolve

sol = solve(system, TRBDF2(linsolve=KLUFactorization()), saveat=1.0)

原页面输出显示 retcode: Success，保存时刻为 0、1、2、3、4、5。这里仅转述文档展示的结果；本任务没有执行代码、计算误差或验证当前依赖组合。

读取各变量的空间解

每个保存时刻的 sol.u[i] 是一个矩阵，行表示变量，列表示网格节点。因此，第三个保存时刻对应 t=2.0；第一行是 Φ，第二行是 Ψ：


sol.u[3]
sol.u[3][1, :]  # 第三个保存时刻所有节点上的 Φ
sol.u[3][2, :]  # 第三个保存时刻所有节点上的 Ψ

文档展示 sol.u[3] 的尺寸为 2×10000 Matrix{Float64}。例如，sol.u[3][1, :] 提取 Φ 在全部 10,000 个节点上的值。变量与矩阵行的对应关系由传入 FVMSystem 的问题顺序确定。

绘制两个变量的时间演化

教程使用 CairoMakie 对每个保存时刻分别绘制 Φ、Ψ 的等值填色图。第一行显示 Φ，第二行显示 Ψ；颜色级别分别取 0:0.1:1 和 1:10:100：


using CairoMakie

fig = Figure(fontsize=38)
for i in eachindex(sol)
    ax1 = Axis(fig[1, i], xlabel=L"x", ylabel=L"y",
        width=400, height=400,
        title=L"\Phi: t = %$(sol.t[i])", titlealign=:left)
    ax2 = Axis(fig[2, i], xlabel=L"x", ylabel=L"y",
        width=400, height=400,
        title=L"\Psi: t = %$(sol.t[i])", titlealign=:left)
    tricontourf!(ax1, tri, sol[i][1, :], levels=0:0.1:1, colormap=:matter)
    tricontourf!(ax2, tri, sol[i][2, :], levels=1:10:100, colormap=:matter)
end
resize_to_layout!(fig)
fig

原文页面包含这组图的输出；本文只保留其文字说明和代码，没有嵌入或生成图片。

无注释版完整示例

以下是原文提供的整段代码，按原结构保留。项目依赖没有在该页面中锁定具体版本：


using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(0, 1, 0, 1, 100, 100, single_boundary=false)
mesh = FVMGeometry(tri)

Φ_bot = (x, y, t, u, p) -> -1 / 4 * exp(-x - t / 2)
Φ_right = (x, y, t, u, p) -> 1 / 4 * exp(-1 - y - t / 2)
Φ_top = (x, y, t, u, p) -> exp(-1 - x - t / 2)
Φ_left = (x, y, t, u, p) -> -1 / 4 * exp(-y - t / 2)
Φ_bc_fncs = (Φ_bot, Φ_right, Φ_top, Φ_left)
Φ_bc_types = (Neumann, Neumann, Dirichlet, Neumann)
Φ_BCs = BoundaryConditions(mesh, Φ_bc_fncs, Φ_bc_types)

Ψ_bot = (x, y, t, u, p) -> exp(x + t / 2)
Ψ_right = (x, y, t, u, p) -> -1 / 4 * exp(1 + y + t / 2)
Ψ_top = (x, y, t, u, p) -> -1 / 4 * exp(1 + x + t / 2)
Ψ_left = (x, y, t, u, p) -> exp(y + t / 2)
Ψ_bc_fncs = (Ψ_bot, Ψ_right, Ψ_top, Ψ_left)
Ψ_bc_types = (Dirichlet, Neumann, Neumann, Dirichlet)
Ψ_BCs = BoundaryConditions(mesh, Ψ_bc_fncs, Ψ_bc_types)

Φ_q = (x, y, t, α, β, γ, p) -> (-α[1] / 4, -β[1] / 4)
Ψ_q = (x, y, t, α, β, γ, p) -> (-α[2] / 4, -β[2] / 4)
Φ_S = (x, y, t, (Φ, Ψ), p) -> Φ^2 * Ψ - 2Φ
Ψ_S = (x, y, t, (Φ, Ψ), p) -> -Φ^2 * Ψ + Φ

Φ_exact = (x, y, t) -> exp(-x - y - t / 2)
Ψ_exact = (x, y, t) -> exp(x + y + t / 2)
Φ₀ = [Φ_exact(x, y, 0) for (x, y) in DelaunayTriangulation.each_point(tri)]
Ψ₀ = [Ψ_exact(x, y, 0) for (x, y) in DelaunayTriangulation.each_point(tri)];

Φ_prob = FVMProblem(mesh, Φ_BCs; flux_function=Φ_q, source_function=Φ_S,
    initial_condition=Φ₀, final_time=5.0)

Ψ_prob = FVMProblem(mesh, Ψ_BCs; flux_function=Ψ_q, source_function=Ψ_S,
    initial_condition=Ψ₀, final_time=5.0)

system = FVMSystem(Φ_prob, Ψ_prob)

using OrdinaryDiffEq, LinearSolve
sol = solve(system, TRBDF2(linsolve=KLUFactorization()), saveat=1.0)
sol.u[3]

sol.u[3][1, :]

using CairoMakie
fig = Figure(fontsize=38)
for i in eachindex(sol)
    ax1 = Axis(fig[1, i], xlabel=L"x", ylabel=L"y",
        width=400, height=400,
        title=L"\Phi: t = %$(sol.t[i])", titlealign=:left)
    ax2 = Axis(fig[2, i], xlabel=L"x", ylabel=L"y",
        width=400, height=400,
        title=L"\Psi: t = %$(sol.t[i])", titlealign=:left)
    tricontourf!(ax1, tri, sol[i][1, :], levels=0:0.1:1, colormap=:matter)
    tricontourf!(ax2, tri, sol[i][2, :], levels=1:10:100, colormap=:matter)
end
resize_to_layout!(fig)
fig