简介一份面向偏微分方程数值模拟的完整实现文档适合具备数学建模和Python编程基础的研究人员、学者及技术爱好者用于流体力学、热传导等复杂物理过程的分析。内容围绕三维立方体域上六个耦合偏微分方程组系统讲解有限差分法的离散化思路、边界条件设定与迭代求解流程并给出完整可运行的Python代码覆盖初始化、主迭代、收敛判断与三维切片可视化。更进一步通过对六个方程的线性组合导出一阶近似方程阐明特定条件下系统可退化为一个扩散方程并针对σ0.1、1、10、100等参数进行数值实验展示不同σ对解形态和扩散过程的影响。资源仅含一个docx文档大小约30KB代码说明与推导过程集中在同一文件内便于查阅和复现。目前已有118人学习浏览适合希望掌握有限差分法并快速上手科学计算和可视化的读者。1. 三维偏微分方程组数值求解为什么有限差分还值得亲手写一遍做论文复现的人迟早会撞上三维偏微分方程组反应-扩散斑图、流体对流、电磁场演化每一类问题都要在三维空间里做数值解法。我见过太多工程师第一反应是导入商业软件或者直接转向 PINNs 求解偏微分方程但当你需要精确控制差分格式、网格密度和精度阶数时最可靠的仍然是亲手写一个有限差分法求解器。这里以 Brusselator 反应-扩散方程组为例完整走一遍基于有限差分法的 Python 实现从一阶近似推导出显式 FTCS 格式到写出可运行的三维求解器和可视化切片。这个方案能解决“复现论文某张三维图但找不到现成代码”的典型诉求也适合把已有二维代码升级到三维的工程师。读完你能直接拿到一份能改、能跑的骨架而不是只能看不能用的片段。2. 一阶近似推导从Taylor展开到三维FTCS差分格式2.1 控制方程选型为什么用Brusselator演示方程组求解在数值解法这个方向上现在能选的路线很多有限差分、有限元、有限体积、谱方法以及这几年火起来的 PINNs。对三维规则网格上的论文复现有限差分法有独特优势实现直观、内存可控、差分公式和论文里的推导一一对应。有限元处理复杂边界更强但三维网格剖分本身会吃掉大量时间PINNs 的优点是免网格但长时间演化精度和可解释性通常比不过传统格式。这里不评价谁更先进只从复制论文结果的角度说有限差分是能在半天内跑出第一张图的方法。三维偏微分方程组的选择会直接影响后面所有代码的复杂度。常见做法是用反应-扩散方程组因为它在斑图动力学论文里出现频率高、非线性项直观、三维效果又足够有辨识度。Brusselator 是其中最常被用来验证数值方法的一个方程写出来是∂u/∂t Du ∇²u a - (b 1)u u²v ∂v/∂t Dv ∇²v bu - u²v其中 u、v 是两个化学组分浓度Du、Dv 是各自扩散系数a、b 是反应速率常数。这里的 ∇² 是三维拉普拉斯算子写作 ∂²u/∂x² ∂²u/∂y² ∂²u/∂z²。选它有三个理由第一它是真正的双变量耦合方程组能体现方程组和单方程的差别第二它有一个简单的均匀稳态解 u0 a、v0 b/a方便在上面叠加扰动来复现论文里的斑图实验第三非线性项 u²v 虽然让隐式求解变得麻烦但显式格式配合足够小的步长可以稳定跑下来正好贴合标题说的“一阶近似”场景。我在实际复现论文时最怕的不是方程难而是找不到参考解。Brusselator 的稳态解、线性稳定性阈值都可以手算哪一步数值出问题立刻能判断是格式写错还是参数选错。如果你要复现的是流体或电磁类 PDEs思路也一样先找一个有解析或半解析结果的简化版本把整套有限差分代码跑稳再替换成目标方程。这也是为什么这篇教程选 Brusselator 而不是一上来就处理 Navier-Stokes 方程组。2.2 时间一阶前向欧拉与空间中心差分的Taylor推导标题里的“一阶近似推导”我的理解分两层时间方向真的只有一阶精度空间方向则是由一阶 Taylor 展开组合出二阶中心差分。先看时间层。将 u(t Δt) 在 t 处做 Taylor 展开u(t Δt) u(t) Δt ∂u/∂t(t) O(Δt²)丢掉 O(Δt²) 项并移项就得到时间导数的前向近似∂u/∂t ≈ (u^{n1} - u^n) / Δt这是最典型的一阶精度格式也叫显式欧拉法。它只依赖当前时刻的 u^n 就能算 u^{n1}不用求解线性方程组代价是时间步长受稳定性限制。空间二阶导数同样从 Taylor 展开出发把 u(x h) 和 u(x - h) 分别展开到二阶项写成u(x h) ≈ u(x) h u(x) h²/2 u(x) u(x - h) ≈ u(x) - h u(x) h²/2 u(x)两式相加后一阶项抵消得到u(x) ≈ (u(x h) - 2u(x) u(x - h)) / h²这就是中心差分。严格说中心差分对空间是二阶精度但推导它的每一步用的都是一阶近似思想——只保留 Taylor 展开的低阶项。很多教程把这种格式整体称作 FTCSForward in Time, Central in Space后面代码里的拉普拉斯算子就是这样写出来的。需要额外说明的是非线性项。比如 u²v在一阶近似框架下直接取当前时刻的 u^n、v^n 代入不做隐式迭代。这种做法相当于把非线性项显式线性化是“一阶近似”在方程层面的体现。它带来的限制是如果 b 很大、反应项时间尺度很短显式格式会被迫用小步长这一点在避坑章节还会遇到。这里要注意一个容易绕晕的地方说“一阶近似”但空间用中心差分是二阶精度那么整体格式算几阶严格说FTCS 是时间一阶、空间二阶整体精度以一阶时间误差为主导所以标题把它归为“一阶近似”并不矛盾。实际论文里如果时间精度不够误差会来自时间方向减小 Δt 能看到改善如果空间精度不够误差不会随 Δt 明显变化。这个特性可以用来快速判断代码里到底是时间步长问题还是网格太粗。2.3 三维显式扩散格式的稳定性限制CFL推导显式 FTCS 格式的稳定性是三维求解器第一个坎。对扩散方程做 von Neumann 稳定性分析时把数值解展开成空间傅里叶模式得到放大因子g 1 2r(cos θx cos θy cos θz - 3)r D Δt / h²要求 |g| ≤ 1 对任意波数成立。观察这个式子括号项最大是 0所有角度都为 0最小是 -6三个角度同时为 π所以最不利的情况给出 1 - 12r ≥ -1即 r ≤ 1/6。换句话说在均匀网格 h dx dy dz 下三维显式扩散格式的时间步长必须满足Δt ≤ h² / (6D)把系数对比一下就能看出维度的影响一维是 1/2二维是 1/4三维是 1/6。维度稳定性限制 r_max说明1D1/2r DΔt/h²2D1/4每个方向各多一个扩散项3D1/6稳定性要求随维度收紧这张表背后是很多人踩过的坑把二维代码的步长直接搬到三维看起来只是加了一个方向实际稳定上限缩水到原来的 2/3。如果你的二维代码刚好用的是临界步长三维必然爆炸。实际操作时一般不会卡在临界值而是取 0.20.5 倍的安全系数。这个系数不是玄学是因为网格越密、反应项越强线性稳定性分析给出的边界就越不够保守。另外需要强调的是上面的分析只针对线性扩散算子。Brusselator 的反应项在 b 较大时会让系统变硬即使满足 r ≤ 1/6反应项也可能把解推成 NaN。所以做论文复现时我的习惯是先关掉反应项单独验证扩散格式再逐步加回非线性项。这样每一步翻车都能快速定位是扩散的问题还是反应的问题。当你把 r ≤ 1/6 写进代码时实际是两步先取 D max(Du, Dv)再用 dt 0.3·h² / (6D)。取最大扩散系数保证两个变量都稳定乘 0.3 而不是 0.5 是因为反应项贡献没有被稳定性分析覆盖。后面第四章主循环的代码就是这么写的。3. 三维网格与差分算子的Python实现向量化是关键3.1 三维Mesh生成与数组索引约定三维有限差分的代码组织第一步是把空间离散化。常见做法是生成三个一维坐标数组再用 meshgrid 展开成三维网格。注意 Python 里数组的维度顺序和空间方向要约定清楚否则后面每一步都在错位。import numpy as np Lx, Ly, Lz 1.0, 1.0, 1.0 Nx Ny Nz 64 # 周期性边界配合均匀网格空间步长用长度除网格数 dx Lx / Nx dy Ly / Ny dz Lz / Nz # indexingij 让数组形状严格等于 (Nx, Ny, Nz) x np.linspace(0.0, Lx, Nx, endpointFalse) y np.linspace(0.0, Ly, Ny, endpointFalse) z np.linspace(0.0, Lz, Nz, endpointFalse) X, Y, Z np.meshgrid(x, y, z, indexingij) print(X.shape) # (64, 64, 64)这段代码的关键在 indexingij。numpy 的 meshgrid 默认用 indexingxy返回的第一个维度是 y 而不是 x在二维画图时差别不大但到三维数组累加、做 np.roll 时会直接导致拉普拉斯算子的 x、y 方向互换。索引约定我统一采用axis0 对应 xaxis1 对应 yaxis2 对应 z数组形状严格等于 (Nx, Ny, Nz)。x、y、z 三个一维数组长度分别是 Nx、Ny、Nz和 linspace 的 endpointFalse 配合保证第 0 个点与最后一个点相邻满足周期边界的需求。如果你熟悉 NumPy可能觉得只建三个一维坐标、不显式创建 X、Y、Z 也能算差分。实际操作中我仍建议创建完整的网格数组因为后面可视化、设置初始扰动、提取切片都会直接用到 X、Y、Z 的维度信息。尤其初始条件如果依赖空间坐标比如 u base 0.1·cos(2πx)cos(2πy)cos(2πz)没有网格数组就得用广播技巧反而更容易出错。网格数组的内存开销在 64³ 量级可以忽略到了 256³ 再优化也不迟。3.2 用np.roll实现周期边界的拉普拉斯算子拉普拉斯算子是三维有限差分代码的核心循环。最容易想到的写法是三层 for 循环遍历所有内部点但那在 Python 里慢到不可接受。常见做法是把差分公式写成数组运算用 np.roll 显式实现“取相邻点”的操作一次调用就完成一个方向的平移。下面是周期边界下的三维拉普拉斯算子def laplacian_periodic(f, dx, dy, dz): # 沿 x 轴的一阶中心差分组合出二阶导数 fxx (np.roll(f, -1, axis0) - 2.0 * f np.roll(f, 1, axis0)) / dx**2 # 沿 y 轴 fyy (np.roll(f, -1, axis1) - 2.0 * f np.roll(f, 1, axis1)) / dy**2 # 沿 z 轴 fzz (np.roll(f, -1, axis2) - 2.0 * f np.roll(f, 1, axis2)) / dz**2 return fxx fyy fzznp.roll(f, -1, axis0) 的意思是沿 x 轴把数组整体左移一格也就是把下标 i1 的元素挪到 i 的位置正好对应中心差分公式里需要 u_{i1} 的地方np.roll(f, 1, axis0) 右移一格对应 u_{i-1}。三个方向各自做一次后相加就等价于三维拉普拉斯。周期性边界自动成立因为 np.roll 会把边界外的元素从头尾绕回来不需要额外处理。这里的性能问题值得单独说。三维数组在 Python 里如果用循环内层每一次都要访问单独元素64³ 网格约 26 万个点三层循环加上反应项计算跑一步要好几秒而 np.roll 版本每步只需几次数组级运算差距在两三个数量级。这也是三维有限差分和二维很大的区别二维还能勉强用循环三维不向量化基本跑不了。如果机器内存紧张可以考虑在原数组上做切片而不是 np.roll但代码可读性会变差我一般先向量化跑通再考虑内存优化。中心差分对光滑解有优势但在三维强对流占主导的问题里会引入振荡届时要考虑迎风差分。这个标题的场景是反应-扩散扩散项刚性为主、对流项缺失中心差分是标准选择。如果你的目标方程组是流体类这里的差分算子和稳定性分析都要重做。3.3 边界条件选择周期边界与Neumann边界的取舍论文复现时边界条件必须跟原论文一致这往往是二维代码升级到三维时最容易被忽略的部分。周期边界实现最简单np.roll 天然支持适合模拟空间均匀、无边界进出的反应-扩散体系但如果论文是封闭体系常见的是零流量 Neumann 边界也就是 ∂u/∂n 0。换边界不只是改边界值差分算子在边界点也需要特殊处理。零流量 Neumann 边界的常见做法是镜像法在计算拉普拉斯算子之前把数组在边界外镜像复制一层。比如对 x 方向令 f[-1, :, :] f[1, :, :] 和 f[0, :, :] f[-2, :, :]然后用同样的中心差分公式计算。等价思路是直接修改差分算子边界节点为单侧差分但镜像法改动最小、不容易引入代码分支。我的建议是先把周期边界跑通验证扩散算子和时间推进都正常再替换成论文要求的边界类型一次只改一个变量。三维数组的边界方向也容易搞混。f[-1, :, :] 是 x 方向的最后一个面f[:, -1, :] 才是 y 方向。每次改边界条件前先验证数组形状和你要处理的轴是否一致。这个问题的坑我在第五章专门记录因为三维里轴写错从报错信息很难看出来往往要等到可视化才发现斑图方向完全反了。代码里的 2.0 * f 会用 float64 计算。三维大规模计算时可以考虑把数组声明成 np.float32内存直接减半但精度也降。是否有必要取决于你是否要算收敛阶如果只是复现斑图案float32 通常够用。4. 主程序时间步进、参数控制与三维可视化4.1 主循环代码与参数表现在把网格、差分算子、控制方程拼成可直接运行的主程序。下面代码从三维网格出发用显式欧拉法推进参数设置参考 Brusselator 论文中常用的量级避免出现需要极小步长的极端参数。import numpy as np # 沿用 3.1 和 3.2 的网格与差分算子 # 参数设置 a, b 1.0, 3.0 Du, Dv 1e-3, 5e-3 # 均匀稳态附近加重幅扰动会翻车控制在小扰动范围 rng np.random.default_rng(42) u np.full((Nx, Ny, Nz), a) 0.05 * rng.uniform(-1.0, 1.0, (Nx, Ny, Nz)) v np.full((Nx, Ny, Nz), b / a) 0.05 * rng.uniform(-1.0, 1.0, (Nx, Ny, Nz)) # 稳定步长 安全系数 0.5 乘三维扩散极限 D_max max(Du, Dv) dt 0.5 * dx**2 / (6.0 * D_max) print(fdt {dt:.6f}, h {dx:.4f}) n_steps 2000 for step in range(n_steps): Lu laplacian_periodic(u, dx, dy, dz) Lv laplacian_periodic(v, dx, dy, dz) # 显式欧拉新值只依赖当前时刻 u_new u dt * (Du * Lu a - (b 1.0) * u u * u * v) v_new v dt * (Dv * Lv b * u - u * u * v) u, v u_new, v_new if step % 500 0: print(fstep {step:5d} u [{u.min():.3f}, {u.max():.3f}] v [{v.min():.3f}, {v.max():.3f}])参数取值不是随便挑的。a1.0、b3.0 让均匀稳态 u01、v03线性稳定性分析给出的阈值在 b 略大于 1a² 时出现因此 b3.0 既有足够的反应驱动又不至于让数值格式太僵。扩散系数 Du 比 Dv 小一个量级这是为了形成经典的 Turing 斑图条件。调试时可以先调大 Du、调小 b降低非线性强度跑通后再逐步逼近论文参数。参数取值作用a1.0反应项常数决定均匀稳态b3.0反应项常数决定振荡/斑图倾向Du1e-3u 的扩散系数Dv5e-3v 的扩散系数h1/64空间步长dt0.5·h²/(6D_max)时间步长含安全系数主循环的逻辑很直接先分别算 u、v 的拉普拉斯算子再计算当前时刻的 u_new、v_new。注意 u_new 的计算用到了 uuv这是 Brusselator 里 u²v 的数组写法。算完更新后用 step % 500 输出一次极值方便确认数值没有发散。初次运行时建议把 n_steps 改成 200看前几步 u 的 min/max 是否保持在小扰动范围。4.2 稳定性判定与步长自动计算三维显式格式的步长不能手拍。我在 2.3 节推出稳定条件后代码里实际写成D_max max(Du, Dv) dt 0.5 * dx**2 / (6 * D_max)第一行取两个扩散系数里的最大值保证对 u、v 都满足稳定条件第二行乘 0.5 留一半安全余量。如果以后改网格只改 Nx、Ny、Nzdt 会自动跟着变小不用每次重新手算。这是三维有限差分里少有的“舒服”之处。如果代码运行中出现了 NaN先别急着检查反应项——大部分情况是时间步长超了稳定性边界。把安全系数从 0.5 降到 0.2 重跑如果 NaN 消失就确认是步长问题如果还在再检查初始扰动是否让 v 变成负数。Brusselator 里 v 出现负值会导致 uuv 变成负值后在反应项里迅速放大这种翻车在参数稍大时非常常见。另一个排查点是边界条件Neumann 边界改错会让边界节点出现异常源项表现是先从边界开始发散。这套参数在 64³ 网格跑 2000 步大多数普通笔记本可以承受。如果网格加到 128³内存和计算量分别涨到 8 倍和 16 倍建议先用 32³ 或 64³ 验证格式确认结果形态正确后再上大网格。三维计算最常见的浪费时间就是网格开太大、步长没配好跑了一半才发现发散重新来。缩小网格是三维调试的后悔药。4.3 切片可视化三维体数据的降维观察三维数组不容易直接可视化。常规做法是取三个中心切片分别看 yz 平面、xz 平面、xy 平面上的浓度分布。这里的数组索引和平面方向要再确认一遍u[:, :, Nz//2] 表示固定 z 轴中面得到的是 x-y 平面u[:, Ny//2, :] 固定 y 中面得到 x-z 平面u[Nx//2, :, :] 固定 x 中面得到 y-z 平面。import matplotlib.pyplot as plt fig, axes plt.subplots(1, 3, figsize(12, 4)) titles [z Lz/2, y Ly/2, x Lx/2] slices [u[:, :, Nz // 2], u[:, Ny // 2, :], u[Nx // 2, :, :]] for ax, data, title in zip(axes, slices, titles): # 转置让数组的第0维对应横轴第1维对应纵轴 im ax.imshow(data.T, originlower, extent[0, 1, 0, 1], cmapviridis) ax.set_title(title) plt.tight_layout() plt.show()imshow 的第一个参数我传了 data.T因为 imshow 把数组第 0 维当作纵轴、第 1 维当作横轴而我们在 3.1 里约定 axis0 是 x传 transpose 后横轴才是 x。originlower 表示 y 轴向上增长extent 把坐标范围设成 0 到 1。这三个参数错一个图里的斑图就会旋转或镜像新手在这里耗掉的时间往往比写主循环还多。切片可视化是快速验证逻辑的手段不是最终成品。要判断结果是否合理可以看初始阶段均匀稳态叠加小扰动后开始的几百步里 u、v 应该保持接近稳态缓慢形成空间结构。如果第一步就出现规则大斑图说明迭代格式或参数有问题。另一个验证点是输出文件如果要做论文对比最好把最后一步的完整三维数组存成 .npy而不是截图后面可以重新切片、算统计量、跑收敛性验证省得重算。n_steps2000 对应物理时间 t_total n_steps·dt。按上面参数估算这个时间尺度足够看到从扰动到初步斑图的过程但离完全稳定可能还不够。论文里常需要更长的演化这时把 n_steps 加大比加大网格更有效因为时间步进是显式的重复利用已经算好的算子即可。5. 避坑记录三维有限差分最容易翻车的4个场景三维有限差分代码写起来比二维多不了几十行但调试难度成倍上升。如果你之前只写过二维有限差分三维的第一感觉会是“代码差不多”第二个感觉才是“哪里都不对”。下面四条是我在实际复现论文和帮同事排查代码时遇到最多的问题按出现频率排序。每一条都按现象、原因、解决的顺序写你遇到类似情况时可以对着查。5.1 内存溢出三维数组比想象中“贵”现象程序运行到一半报 MemoryError或者系统开始疯狂使用交换分区计算速度掉到原来十分之一。原因三维数组的元素数量是三个维度相乘64³ 就是 262144 个点两个变量用 float64 也就 4 MB本身不贵。但 np.roll 每次都会创建临时数组拉普拉斯算子一次要生成 6 份临时数组2000 步循环里的临时对象加在一起即使被垃圾回收峰值内存也会翻好几倍。如果网格加到 128³临时数组峰值轻松突破几百 MB内存不够就翻车。解决先降网格到 32³ 或 64³其次可以把三个方向的差分合并成一次数组切片运算减少临时数组数量再不够就把数组类型改成 np.float32。我一般优先改网格尺寸因为调试期不需要高分辨率小网格能把一轮试验时间从小时级压到分钟级。5.2 时间步长失控导致的数值爆炸现象连续几步后 u、v 的 max 值从几跳到 1e20随后全是 nan或者极值在某个迭代步突然变成 -nan。原因九成是时间步长超过三维稳定极限。很多人从二维代码改三维时只改了拉普拉斯算子的维度忘了稳定条件已经收紧。剩下的一成是初始扰动太大v 在扰动后被推成负值反应项 u²v 瞬间产生巨大负反馈显式格式救不回来。解决先按 dt 0.2·dx² / (6·D_max) 重跑确认爆炸消失后逐步提高安全系数再把初始扰动幅度从 0.05 降到 0.01这一步可以排除反应项问题。如果两个都调完仍然爆炸才考虑是不是方程系数写错。这里没有太多玄学绝大多数数值爆炸都是步长问题。5.3 索引不一致三维数组放反坐标轴现象数值上没报错但可视化时看到的斑图方向跟预期完全相反或者在 x 方向看到本应在 y 方向的条带。原因meshgrid 默认是 indexingxy如果你图省事没传 indexingij数组的第 0 维就变成了 y。后续所有代码都按 axis0 是 x 来写等于把 x、y 整体换位。拉普拉斯算子在均匀网格下还能算但初始条件和边界条件的方向会全都错位。这个 bug 最隐蔽的点在于均匀网格下扩散方程对方向不敏感算出来的解本身没问题只有等你加各向异性项或者对照论文对比时才暴露。解决从一开始统一用 indexingij并且在文件开头写一行注释固定 axis0x、axis1y、axis2z。每次加新代码前先打印一个 2×2×2 的小数组检查维度顺序。三维数组的轴顺序就像方向盘的左右养成固定习惯后基本不会再犯。5.4 反应项与扩散项刚性冲突现象满足扩散稳定性条件后结果仍然剧烈振荡或时间步长被迫取得极小2000 步耗时比预期多一个量级。原因Brusselator 的线性稳定性分析显示b 超过阈值后均匀态失稳当 b 进一步增大反应项的特征时间远小于扩散项显式格式要同时满足两者的稳定限制实际上被反应项拖着走。扩大 b 想让斑图更明显结果步长越来越小这就是典型的刚性冲突。解决先确认参数在论文给定的范围内不要随手加大 b其次把安全系数降到 0.2如果步长实在太小可以只对反应项做隐式处理扩散项保持显式也就是算子分裂。u 的更新可以写成一行u_new (u dt * (Du * lap(u) a u*u*v)) / (1 dt * (b 1.0))v 的方程思路一样。反应项是局部的每个网格点独立隐式处理不需要解大规模线性方程组比全隐式简单得多。我一般先用小 b 验证代码逻辑再逼近目标参数避免一开始就被刚性问题带偏。以上四条有一个共同点问题都不是出在公式推导上而是出在三维对维度、内存、步长的放大效应上。二维代码能容忍的不规范在三维会被放大成一个小时的无效等待。所以我会在每次换维度或换网格时先把这些检查项过一遍再开始跑长任务。6. 进阶用解析解验证收敛阶把一阶时间格式换成RK26.1 验证技巧三维扩散方程解析解与L2误差计算先验证再改格式是避免重复劳动的关键。去掉反应项后三维扩散方程在周期边界下有解析解取初值为 cos(kx)cos(ky)cos(kz)解随时间按因子 exp(-3Dk²t) 衰减。固定 dt、h 跑若干步后与解析解比较计算 L2 误差把网格从 16³ 加密到 32³、64³误差应接近按二阶下降。这个验证只需要改一个参数、加几行计算能一次性确认网格生成、拉普拉斯算子、边界条件三个环节都正确。我在换机器或换 NumPy 版本后会先跑这个用例当作回归测试。6.2 把前向欧拉升级为RK2的改法扩散时间格式从欧拉换成二阶 Runge-Kutta 很简单算子完全复用只在步进函数上做改动。lap lambda F: laplacian_periodic(F, dx, dy, dz) for step in range(n_steps): # 中间步半步欧拉 u_mid u 0.5 * dt * (Du * lap(u) a - (b 1.0) * u u * u * v) v_mid v 0.5 * dt * (Dv * lap(v) b * u - u * u * v) # 完整步用中间时刻的斜率 u_new u dt * (Du * lap(u_mid) a - (b 1.0) * u_mid u_mid * u_mid * v_mid) v_new v dt * (Dv * lap(v_mid) b * u_mid - u_mid * u_mid * v_mid) u, v u_new, v_new逻辑说明先做半步欧拉得到中点值 u_mid、v_mid再用中点斜率做完整步更新。相比前向欧拉RK2 对反应项的相位误差更小稳定性也有改善但稳定步长不建议因此放大保守仍是好习惯。这里的 lap 就是第三章写的拉普拉斯算子。我在复现论文时一般先保留欧拉格式跑出定性图案确认参数区间和边界条件都对了再换 RK2 提高时间精度配合 6.1 的解析解验证把误差压下去。最后补一点教训三维有限差分的精度验证一定在换大网格之前做否则错误会被高分辨率浪费掉的算力掩盖。希望帮到你。本文还有配套的精品资源点击获取