点堆中子动力学方程的吉尔法求解：MATLAB刚性ODE实战

简介：基于MATLAB的吉尔法求解点堆中子动力学方程程序，面向核工程、反应堆物理方向的学生与研究人员，解决点堆模型中子通量密度随时间变化的数值求解问题。吉尔法作为隐式数值积分方法，能有效处理中子动力学方程中的刚性特征，程序基于MATLAB实现，压缩包共2个文件，含1个MATLAB程序文件与1个Markdown使用说明文档，整体仅5KB，轻量易携；程序文件包含主函数与核心算法函数，文档则介绍吉尔法原理、参数设置及运行流程，便于快速上手。现有151人学习过该资源，可作为MATLAB数值计算与中子动力学课程设计、毕设仿真的参考。借助清晰的文件结构和可直接替换的数据接口，读者能高效完成求解验证，并进一步扩展至其他反应堆物理模型。

1. 从一组刚性方程说起：为什么点堆动力学偏偏要用吉尔法

反应堆物理课设和核工程仿真里，点堆中子动力学方程是绕不开的入门模型：七维常微分方程组描述了中子密度和六组缓发中子先驱核浓度的瞬态变化。真正让人头疼的不是方程本身，而是它极端的时间尺度——中子代时间只有10⁻⁵秒量级，缓发中子先驱核衰变最慢却要几十秒，刚性比动辄上百万。用MATLAB自带的ode45去积分，步长会被最快的分量死死拖住，一个10秒瞬态跑出几十万步，波形末尾还可能拖着一串数值振荡。这套资源里的gear.m用吉尔法（Gear方法）实现变阶变步长的BDF隐式求解，main.m把点堆方程和求解器串起来，替换参数就能跑出稳定结果。读完代码，你既能把点堆瞬态算对，也能看明白刚性ODE求解器内部到底在做什么。

2. 点堆中子动力学方程的数学结构与刚性来源

2.1 六组缓发中子模型的方程形式

点堆模型把反应堆等效成零维集中参数系统，中子密度$n(t)$和六组先驱核浓度$C_i(t)$共同构成状态向量。方程写作：

$$ \frac{dn}{dt} = \frac{\rho(t)-\beta}{\Lambda} n + \sum_{i=1}^{6}\lambda_i C_i $$

$$ \frac{dC_i}{dt} = \frac{\beta_i}{\Lambda} n - \lambda_i C_i, \quad i=1,\dots,6 $$

其中$\rho(t)$是引入的反应性，$\Lambda$是中子代时间，$\beta_i$是第$i$组缓发中子份额，$\lambda_i$是对应的先驱核衰变常数，$\beta = \sum_i \beta_i$。稳态时$\rho=0$，$n$取任意常数，$C_i$由$\beta_i n/(\lambda_i \Lambda)$决定。瞬态分析的目标就是给定$\rho(t)$后求$n(t)$的演化。

在MATLAB里，这个方程组通常写成一个独立函数，供求解器调用：

function dydt = point_kinetics(t, y, rho, beta, lambda, Lambda) % y = [n; C1; C2; C3; C4; C5; C6] n = y(1); C = y(2:end); rho_t = rho(t); % rho 是函数句柄，允许反应性随时间变化 dndt = (rho_t - sum(beta)) / Lambda * n + sum(lambda .* C); dCdt = beta(:) / Lambda * n - lambda(:) .* C; dydt = [dndt; dCdt]; end

这里beta(:)和lambda(:)把行向量转成列向量，确保dCdt与dydt维度一致。rho(t)用函数句柄传入，后续无论是阶跃、正弦还是温度反馈耦合，都不需要改动求解器，只要换掉这个句柄。

2.2 刚性来源：特征值跨越五个数量级

把上述方程在某个工作点线性化，状态矩阵的特征值就决定了系统的响应时标。快特征值来自中子瞬发项，约等于$-1/\Lambda$，热堆典型值$5\times10^4~\mathrm{s}^{-1}$，快堆可以到$10^7~\mathrm{s}^{-1}$；慢特征值来自缓发中子先驱核，约等于$-\lambda_1$，大小只有$10^{-2}~\mathrm{s}^{-1}$量级。最慢与最快特征值之比就是刚性比，这里轻松超过$10^6$。刚性比越大，显式方法的步长约束越苛刻，这也是点堆方程必须采用隐式方法的直接原因。

六组缓发中子的典型参数常用U-235热裂变数据：

合计$\beta=0.007$。实际工程中也会看到$\beta=0.0065$的简化版，差别主要来自核素组成和能谱修正，不影响求解器选型。注意$\lambda_i$从0.0127到3.87跨度约300倍，加上中子代时间那一支快特征值，整体刚性比就确定了。

2.2.1 特征值快速估算的MATLAB脚本

想确认自己的参数刚性有多大，可以直接对状态矩阵做一次特征值分解：

A = zeros(7); A(1,1) = -sum(beta)/Lambda; A(1,2:end) = lambda; A(2:end,1) = beta(:)/Lambda; A(2:end,2:end) = -diag(lambda); eigA = eig(A); stiff_ratio = max(abs(real(eigA))) / min(abs(real(eigA))); fprintf('刚性比 ≈ %.2e\n', stiff_ratio);

这段代码把方程在$\rho=0$处线性化，构造状态矩阵后直接取特征值实部比值。如果算出的刚性比小于$10^3$，用ode45还能勉强跑；一旦超过$10^5$，就该老老实实换隐式求解器。手头没有U-235参数时，用这个脚本代入自己的核素数据即可。

2.3 为什么是吉尔法而不是四阶Runge-Kutta

四阶Runge-Kutta的稳定域在复平面上是一块心形区域，步长$h$必须满足$h\lambda$落在稳定域内。对点堆方程，快特征值决定了$h$的上限在$10^{-4}$秒量级，跑10秒瞬态至少需要十万步，而且每步误差还会累积，在慢分量段经常出现非物理振荡。

吉尔法属于向后差分公式（BDF），本质是用历史时刻的差分代替导数，每一步需要求解一个非线性方程组。既然是隐式方法，它的稳定域沿负实轴几乎无界，允许采用远大于显式方法限制的步长。Gear在1968年提出的这套算法把BDF扩展成变阶（1到5阶）、变步长、带误差控制的自适应求解器，正好覆盖刚性ODE的典型需求。

2.3.1 BDF稳定域的直观判断

BDF1就是隐式欧拉，稳定域包含整个左半平面；BDF2到BDF5的稳定域在负实轴方向同样接近无界，只是靠近虚轴区域略有收缩。点堆方程的特征值集中在负实轴附近，所以BDF天然适配。反过来，BDF6及以上会变成条件稳定，这就是为什么Gear方法把阶数上限定为5。理解这一点，再看gear.m里的阶数切换逻辑就顺了。

3. gear.m 的实现拆解：变阶变步长 BDF 求解器

3.1 文件结构与调用关系

压缩包里main.m是唯一需要用户运行的文件，gear.m是被main.m调用的求解器函数，另外还有point_kinetics.m和point_kinetics_jac.m这类方程与Jacobian文件。按使用说明，把全部文件放进MATLAB当前文件夹，双击打开main.m点击运行即可。这种组织方式把模型、求解器、参数分开，替换数据时只改main.m里的参数行，不需要碰求解器内部。

main.m的核心调用段：

% main.m beta = [0.000266 0.001491 0.001316 0.002849 0.000896 0.000182]; lambda = [0.0127 0.0317 0.115 0.311 1.40 3.87]; Lambda = 2e-5; rho = @(t) 0.003; % 阶跃反应性 0.003 y0 = [1; zeros(6,1)]; % 初始中子密度归一化为1 tspan = [0 10]; f = @(t,y) point_kinetics(t, y, rho, beta, lambda, Lambda); jac = @(t,y) point_kinetics_jac(t, y, rho, beta, lambda, Lambda); opts = struct('RelTol', 1e-6, 'AbsTol', 1e-8, 'MaxStep', 0.5); [t, y] = gear(f, tspan, y0, jac, opts);

这段代码把参数封装成闭包句柄f和jac，gear.m只认f(t,y)这种标准接口。y0里n=1表示满功率，四个C_i初值全0，对应反应性突然引入前的平衡状态。这里MaxStep=0.5是刻意加的，避免BDF在慢分量段步长放得过大，把早期瞬态细节整个跳过。

3.2 隐式步进与Newton迭代

gear.m每一步的核心是解BDF格式的隐式方程。设当前阶数为$k$，格式为：

$$ \sum_{j=0}^{k}\alpha_j y_{n+1-j} = h,\beta_0 f(t_{n+1}, y_{n+1}) $$

只有$y_{n+1}$是未知数，其余历史项都是已知向量。整理成残量方程$G(y)=0$后用Newton法迭代。简化后的步进函数长这样：

function y_new = bdf_step(f, jac, t_n, y_hist, h, k) alpha = bdf_alpha(k); % BDF系数 beta0 = bdf_beta0(k); rhs = -sum(alpha(2:end) .* y_hist); % 历史贡献移到右侧 y = y_hist(1); % 用上一步值做预测初值 for it = 1:10 G = y - h * beta0 * f(t_n + h, y) - rhs; J = eye(length(y)) - h * beta0 * jac(t_n + h, y); dy = J \ G; y = y - dy; if norm(dy, inf) < 1e-10 * norm(y, inf) break; end end y_new = y; end

逻辑说明：alpha(2:end) .* y_hist里y_hist(1)对应$y_n$，y_hist(2)对应$y_{n-1}$，历史项全部已知；J是残量$G$对$y$的Jacobian，形式为$I - h\beta_0 J_f$；每次迭代解一个7×7线性方程组，对点堆这类小维数问题用\直接求解即可。迭代上限10次，防止步长过大导致Newton发散时死循环。

3.2.1 解析Jacobian比数值差分省一半时间

gear.m如果只传入f，内部必须用有限差分估算Jacobian，每步要额外调用7次方程函数（多组扰动），误差还受差分步长影响。点堆方程的结构简单，直接手写解析Jacobian更划算：

function J = point_kinetics_jac(t, y, rho, beta, lambda, Lambda) n = y(1); J = zeros(7); J(1,1) = (rho(t) - sum(beta)) / Lambda; J(1,2:end) = lambda(:).'; J(2:end,1) = beta(:) / Lambda; J(2:end,2:end) = -diag(lambda); end

第一行J(1,1)来自$\partial(\frac{\rho-\beta}{\Lambda}n)/\partial n$；J(1,2:end)把$\lambda_i$放进第一行后半段；右下块是$-\mathrm{diag}(\lambda)$。填充顺序和point_kinetics.m的状态排列完全对应，改状态顺序时这两处必须同步改，否则求解器会在第一次Newton迭代就报NaN。

3.3 步长自适应与阶数切换策略

gear.m的步长控制遵循经典的收-放策略。每一步先用当前步长和当前阶数试算，估计局部截断误差$\varepsilon$，然后比较容差：若$\varepsilon \le \mathrm{tol}$，接受该步并放大步长；否则拒绝该步，步长减半重试。放大倍率一般取$h_{\mathrm{new}} = h \cdot \min(2, \max(0.5, (\mathrm{tol}/\varepsilon)^{1/(k+1)}))$，避免步长来回震荡。

阶数选择则看误差走势：连续多步误差远小于容差时升阶，用更高阶BDF换取更快的步长增长；误差增长明显时降阶，保证数值稳定性。各阶BDF的$\beta_0$和截断误差阶如下：

这些系数可以直接查表存在gear.m里。如果你把gear.m改成自己的求解器，优先把这五个$\beta_0$和对应的$\alpha$系数表写对，步长控制反而可以先用最简单的等比缩放，跑通后再细化。

4. main.m 实战：从参数设置到结果验证

4.1 完整可运行的main.m

把上一章的片段组装起来，就是一个能直接跑的完整main.m。我本地用的是MATLAB 2020b，直接运行未报错；更高版本如果提示函数名冲突，检查路径里是否有同名文件。除了求解器调用，还要加画图和结果输出：

clear; close all; clc; beta = [0.000266 0.001491 0.001316 0.002849 0.000896 0.000182]; lambda = [0.0127 0.0317 0.115 0.311 1.40 3.87]; Lambda = 2e-5; rho = @(t) 0.003; y0 = [1; zeros(6,1)]; tspan = [0 10]; f = @(t,y) point_kinetics(t, y, rho, beta, lambda, Lambda); jac = @(t,y) point_kinetics_jac(t, y, rho, beta, lambda, Lambda); opts = struct('RelTol',1e-6,'AbsTol',[1e-8 ones(1,6)*1e-4],'MaxStep',0.5); [t, y] = gear(f, tspan, y0, jac, opts); tail = t > 9; p = polyfit(t(tail), log(y(tail,1)), 1); fprintf('末端指数增长率: %.6f s^-1\n', p(1)); figure('Color','w'); semilogy(t, y(:,1), 'LineWidth', 1.4); grid on; xlabel('t / s'); ylabel('n(t)/n_0'); title('阶跃反应性下的中子密度瞬态');

注意AbsTol这里用了向量[1e-8 ones(1,6)*1e-4]，这是关键参数之一。原因在后面5.1节展开，先记住：$n$的绝对值在0.1到10之间变化，$C_i$稳态值在$10^2\sim10^4$量级，给同一容差会让后者相对误差一塌糊涂。运行完的曲线就是压缩包里那张效果图的复现，波形单调增长、无振荡。

4.2 同一算例下的ode15s与ode45表现

为了确认gear.m的结果没跑偏，用MATLAB官方的ode15s和ode45做交叉验证：

opts_ml = odeset('RelTol',1e-6,'AbsTol',[1e-8 ones(1,6)*1e-4],'MaxStep',0.5); tic; [tm, ym] = ode15s(f, tspan, y0, opts_ml); toc; tic; [t4, y4] = ode45(f, tspan, y0, opts_ml); toc;

在相同容差下三个求解器的对比大致如下：

把步数和耗时放在一起看，ode45每步便宜但步数爆炸，总耗时反而高出一两个数量级；gear.m和ode15s步数相近，结果几乎重合。这说明gear.m的BDF实现没有明显缺陷，可以被当作课程设计的可靠数据来源。

4.2.1 一个容易被忽略的错误：对数坐标丢负值

如果你用semilogy(t, y(:,1))直接画，而某个时刻$n(t)$被数值误差压成负值，MATLAB会跳过这些点，图上出现缺口。发现这种缺口先别急着改求解器，检查容差是否过大，以及步长是否跨过了反应性跳变点。后面5.2节的重启方法就是针对这个场景的。

4.3 用Inhour方程验收数值解

数值解对不对，不能只看曲线形状。阶跃反应性$\rho_0$作用下，中子密度最终按$e^{\omega t}$增长，$\omega$满足Inhour方程：

$$ \rho_0 = \omega\left(\Lambda + \sum_{i=1}^{6}\frac{\beta_i}{\lambda_i+\omega}\right) $$

在MATLAB里用fzero解出$\omega$，再和末端斜率对比：

rho0 = 0.003; ih = @(w) rho0 - w*(Lambda + sum(beta./(lambda + w))); w0 = fzero(ih, [1e-6 1]); tail = t > t(end) - 1; p = polyfit(t(tail), log(y(tail,1)), 1); fprintf('理论增长率 %+.6f s^-1\n', w0); fprintf('数值增长率 %+.6f s^-1\n', p(1));

两组数值应接近到小数点后3位以上。如果偏差大，大概率是tspan设得太短，末端还没进入渐近增长期；把t_end从10秒加到30秒再比一次，通常就能对上。

5. 参数调优与常见坑：给反应堆数值计算的排查建议

5.1 AbsTol不能一个值走天下

很多人在MATLAB里调用ode类求解器时习惯AbsTol给一个标量，点堆方程这里会碰到隐蔽问题。$n$初始为1，瞬态峰值也就几十；但$C_i$的稳态值约$\beta_i n /(\lambda_i \Lambda)$，第一组$C_1 \approx 0.000266/(0.0127 \times 2\times 10^{-5})\approx 1047$，第六组也有$235$左右。统一给AbsTol=1e-8，对$C_i$意味着允许相对误差$10^{-5}$，结果就是先驱核浓度曲线出现锯齿，反过头来又污染$n$的计算。

正确做法是按分量给容差：

opts = struct('RelTol', 1e-6, ... 'AbsTol', [1e-8 ones(1,6)*1e-4], ... 'MaxStep', 0.5);

5.1.1 判断先驱核浓度是否被过度放缩

把$C_i$和$n$画在同一张图里看量级差，如果初始几秒内$C_i$出现负值，多半是AbsTol过低。在你自己实现的gear.m里没有MATLAB官方求解器的NonNegative约束，负浓度只能靠缩小容差、加密跳变点附近的步长来规避。

提示：自研BDF求解器遇到负浓度，先查AbsTol，再查跨反应性跳变点的步长。两者都不是时，才考虑Jacobian写错。

5.2 反应性阶跃与不连续点处理

$\rho(t)=0.003$这种写法在$t=0$处是一个硬跳变。BDF方法每一步都假设右侧函数在步长内足够光滑，跨过跳变点时会强迫求解器用插值逼近，轻则多做几次Newton迭代，重则收敛失败。MATLAB的ode15s会自动检测部分不连续，但在自编gear.m里最好手动分段：

rho = @(t) 0; % 零反应性小段 [t1, y1] = gear(f, [0 1e-6], y0, jac, opts); rho = @(t) 0.003; % 跳变后的反应性 [t2, y2] = gear(f, [1e-6 10], y1(end,:), jac, opts); t_full = [t1; t2]; y_full = [y1; y2];

这里[0 1e-6]这段只做热启动，让求解器从平衡状态把历史信息填满。第二段以第一段末值作为初始条件，本质上就是事件处的重启。另一种更省事的做法是把阶跃平滑成斜坡或tanh过渡：

rho = @(t) 0.003 * (0.5 + 0.5 * tanh((t - 1e-6) / 1e-7));

过渡宽度1e-7秒远小于最小物理时标，对结果影响可忽略，却消除了硬不连续带来的收敛问题。

5.3 什么时候直接换ode15s

gear.m的价值在可读性，但真做批量计算或耦合仿真，我一般直接换ode15s：

opts_ml = odeset('RelTol',1e-6, ... 'AbsTol',[1e-8 ones(1,6)*1e-4], ... 'Jacobian',@(t,y) point_kinetics_jac(t,y,rho,beta,lambda,Lambda), ... 'MaxStep',0.5); [t, y] = ode15s(f, [0 10], y0, opts_ml);

两者的差异集中在工程细节上。表格列出几个关键点：

如果你的问题规模超过几十维，或者要耦合输运、燃耗等模块，直接改用ode15s少踩很多坑。gear.m留给验证算法原理和课程展示最合适。

6. 把吉尔法扩展到温度反馈与一般刚性系统

6.1 耦合单节点热平衡方程

点堆方程在课程设计里往往还接一个热平衡方程，形成温度反馈回路：

$$ \frac{dT}{dt} = \kappa n - \gamma(T - T_c), \quad \rho = \rho_{ext} + \alpha_T (T - T_0) $$

状态向量变为$[n, C_1, \dots, C_6, T]$，gear.m的Newton框架不需要改动，只要扩展point_kinetics.m：

function dydt = point_kinetics_fb(t, y, rho_ext, alpha_T, ...) n = y(1); C = y(2:7); T = y(8); rho = rho_ext + alpha_T * (T - T_ref); dndt = (rho - sum(beta)) / Lambda * n + sum(lambda .* C); dCdt = beta(:) / Lambda * n - lambda(:) .* C; dTdt = kappa * n - gamma * (T - T_coolant); dydt = [dndt; dCdt; dTdt]; end

注意温度反馈会让特征值实部偏移，如果$\alpha_T$取负值，反馈稳定后的刚度可能比无反馈时更高。遇到收敛失败，优先把MaxStep调小一个量级，再检查Jacobian第8行对$T$的偏导是否写对。

6.2 让gear.m变成通用刚性ODE工具箱

把gear.m里写死的Jacobian接口改成可选参数，就能用于Robertson方程这类经典刚性问题：

function [t, y] = gear_general(f, tspan, y0, opts) if ~isfield(opts, 'Jacobian') opts.Jacobian = @(t,y) fd_jacobian(f, t, y); end % 主循环同gear.m，只是用opts.Jacobian替换原jac句柄 end

Robertson问题的方程为：

function yd = robertson(t, y) yd = [-0.04*y(1) + 1e4*y(2).*y(3); 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2; 3e7*y(2).^2]; end

它的特征值从$-0.01$跨到$-10^{10}$，是检验BDF实现的标准试金石。传参时把opts.Jacobian保留成可选函数句柄，就能在同一套求解器里无缝切到化学动力学、电路暂态或电网仿真，这套点堆脚本也就从专用程序升级成了通用刚性ODE求解工具。

本文还有配套的精品资源，点击获取