在 Julia 中模拟 Keller–Segel 趋化模型并制作双场动画-未完纪

本文译写自 SciML 的 FiniteVolumeMethod.jl 官方教程《Keller-Segel Model of Chemotaxis》，并对代码中的通量、数组形状和动画联动补充说明。原页未署个人作者名，归属为 FiniteVolumeMethod.jl 文档维护者；页面由 Literate.jl 生成。

核验时读取的是 stable 文档。页面设置区写明：使用 Documenter.jl 1.8.0、Julia 1.11.2，于 2025 年 1 月 1 日生成。这是文档构建日期，不是首次发表日期；页面也没有锁定 FiniteVolumeMethod.jl 及全部依赖的具体版本。以下源码与显示输出均来自原文，本次只做静态审查，没有运行求解或录制动画。

两个变量共同组成趋化系统

本例考虑的 Keller–Segel 模型为：

∂u/∂t = ∇²u − ∇·[(c u / (1 + u²)) ∇v] + u(1 − u)
∂v/∂t = D ∇²v + u − a v

区域是正方形 [0,100]²，两个变量采用齐次 Neumann 边界条件。需要分别为 u 和 v 定义问题，然后组合成一个系统。

编者导读：把 u 理解为细胞密度、v 理解为化学信号时，第一式同时包含扩散、受信号梯度驱动的趋化和 logistic 增长；第二式包含扩散、由 u 产生的信号以及线性衰减。下面的实现明确保留了两个变量之间的耦合，没有把它们当成互不相关的两个扩散方程。

原创技术示意图：上半部分展示通量与源项的关系，下半部分展示同一 Observable 如何同步更新双场图。不是模拟结果或动画截图。

网格、边界与通量

下面保留原文建立问题的完整代码。先在正方形上建立 250×250 的节点网格并转换成 FVMGeometry，再为两个变量分别建立边界条件和通量函数：

using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(0, 100, 0, 100, 250, 250, single_boundary=true)
mesh = FVMGeometry(tri)
bc_u = (x, y, t, (u, v), p) -> zero(u)
bc_v = (x, y, t, (u, v), p) -> zero(v)
BCs_u = BoundaryConditions(mesh, bc_u, Neumann)
BCs_v = BoundaryConditions(mesh, bc_v, Neumann)
q_u = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    u = αu * x + βu * y + γu
    ∇u = (αu, βu)
    ∇v = (αv, βv)
    χu = p.c * u / (1 + u^2)
    _q = χu .* ∇v .- ∇u
    return _q
end
q_v = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    ∇v = (αv, βv)
    _q = -p.D .* ∇v
    return _q
end
S_u = (x, y, t, (u, v), p) -> begin
    return u * (1 - u)
end
S_v = (x, y, t, (u, v), p) -> begin
    return u - p.a * v
end
q_u_parameters = (c=4.0,)
q_v_parameters = (D=1.0,)
S_v_parameters = (a=0.1,)
u_initial_condition = 0.01rand(DelaunayTriangulation.num_points(tri))
v_initial_condition = zeros(DelaunayTriangulation.num_points(tri))
final_time = 1000.0
u_prob = FVMProblem(mesh, BCs_u;
    flux_function=q_u, flux_parameters=q_u_parameters,
    source_function=S_u,
    initial_condition=u_initial_condition, final_time=final_time)
v_prob = FVMProblem(mesh, BCs_v;
    flux_function=q_v, flux_parameters=q_v_parameters,
    source_function=S_v, source_parameters=S_v_parameters,
    initial_condition=v_initial_condition, final_time=final_time)
prob = FVMSystem(u_prob, v_prob)

原文创建系统时显示：

FVMSystem with 2 equations and time span (0.0, 1000.0)

这段代码的关键是 q_u 接收两个变量的局部线性表示系数。它用 u = αu*x + βu*y + γu 重建当前位置的密度，得到自身梯度 ∇u = (αu,βu) 和信号梯度 ∇v = (αv,βv)，再计算状态依赖的趋化系数 χu = c*u/(1+u²)。

代码返回的密度通量为 q_u = χu∇v − ∇u；信号通量则是 q_v = −D∇v。按照“时间导数 = −通量散度 + 源项”的写法，就得到开头的两条方程。密度源项 u*(1-u) 给出 logistic 增长，信号源项 u-a*v 同时表示产生与衰减。

原文参数为 c=4.0、D=1.0、a=0.1。密度初值使用 0.01rand(...) 生成小幅随机扰动，信号初值全为零，终止时间是 1000.0。两个 FVMProblem 共用网格，并按 FVMSystem(u_prob, v_prob) 的顺序组成双变量系统。

求解并同步显示两个场

模型定义完毕后，教程使用 Sundials 的 CVODE_BDF，将线性求解器设为 GMRES，每隔一个时间单位保存一次结果。然后用 CairoMakie 建立并排的两个热图，把它们绑定到同一个帧索引：

using OrdinaryDiffEq, Sundials, CairoMakie
sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), saveat=1.0, parallel=Val(false)) 
fig = Figure(fontsize=44)
x = LinRange(0, 100, 250)
y = LinRange(0, 100, 250)
i = Observable(1)
axu = Axis(fig[1, 1], width=600, height=600,
    title=map(i -> L"u(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
axv = Axis(fig[1, 2], width=600, height=600,
    title=map(i -> L"v(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
u = map(i -> reshape(sol.u[i][1, :], 250, 250), i)
v = map(i -> reshape(sol.u[i][2, :], 250, 250), i)
heatmap!(axu, x, y, u, colorrange=(0.0, 2.5), colormap=:turbo)
heatmap!(axv, x, y, v, colorrange=(0.0, 10.0), colormap=:turbo)
resize_to_layout!(fig)
record(fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol);
    framerate=60) do _i
    i[] = _i
end;

这里 i = Observable(1) 是两张图共用的控制量。每次改变 i[]，左右标题会从同一个 sol.t[i] 取时间，左图使用 sol.u[i][1,:]，右图使用 sol.u[i][2,:]。它们分别被重塑为 250×250 数组，与前面的节点网格相对应。

u 的颜色范围固定为 (0.0,2.5)，v 的颜色范围固定为 (0.0,10.0)，两者都用 :turbo 色图。固定范围使每一场在不同时间的颜色有一致尺度，但两张图的范围不同，不能把相同颜色直接理解成相同数值。

record() 遍历保存的解索引，每一步令 i[] = _i，以每秒 60 帧录制 MP4。saveat=1.0 决定保存的模拟时间间隔，framerate=60 决定视频播放速率，它们表达的是不同的时间尺度。

原文在录制代码之后提供了一段动画，并感叹生成的图案很精彩。可到原教程页面查看。本文没有下载或重新生成该动画，不把示意图当作实际数值结果。

运行前需要注意的边界

以下是本次静态审核补充，不是原文的运行验证结果：

依赖版本：教程包含 FiniteVolumeMethod、DelaunayTriangulation、OrdinaryDiffEq、Sundials 和 CairoMakie，但没有提供完整锁定环境。文档当时使用的 Julia 版本不能替代包版本约束，也不能保证现行包组合直接兼容。
随机性：源文没有固定随机种子，因此重复运行不应预期得到逐像素相同的图案。本文保留原始初始化方式，没有自行添加种子后声称结果已复现。
资源需求：250×250 节点对应 62,500 个位置、两个状态分量。较长时间范围和逐时保存会占用明显的内存与计算资源，动画渲染还会增加开销。本次没有测量耗时、峰值内存或收敛表现。
形状约定：reshape(...,250,250) 和坐标 LinRange(...,250) 与原网格相配。如果修改节点数或网格形式，不能仍照抄这些尺寸或假定节点顺序不变。
输出文件：原路径 joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4") 以脚本目录为基准，适用于原教程目录结构。运行前需要确认目标目录存在、可写，并检查同名文件，避免覆盖需要保留的视频。本文没有创建目录或写入视频。
求解结果：源码未展示求解返回码、误差或网格收敛性检查；不能仅凭能绘图或原文的系统摘要，就断言解已经满足某个精度标准。