Runge-Kutta 法

一、直观理解：一步之内多点取样再加权

Runge–Kutta（RK）方法的想法是：在一个步长 $h$ 内部多处“试探”斜率，再把这些斜率加权平均，得到更接近真实平均变化率的前进步长。

经典四阶 RK（RK4）常用四个斜率：

$k_1$ ：起点斜率；
$k_2$ ：按 $k_1$ 走到区间中点附近的斜率；
$k_3$ ：按 $k_2$ 再估一次中点斜率；
$k_4$ ：按 $k_3$ 估终点斜率；

最终

y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4).

这像测量一条弯曲河道的平均流向：不是只在码头测一次，而是在河中多个位置测流速再综合。因此在相同步长下，RK4 通常远比 Euler、改进 Euler 准确。

MATLAB 中最常用的 ode45 属于 Dormand–Prince 自适应 RK：同时计算 4 阶与 5 阶公式，用其差估计局部误差并自动调整 $h$ 。建模竞赛里，非刚性 ODE 的默认首选几乎总是 ode45。

二、适用场景与优缺点

类型	说明
适用场景	非刚性或弱刚性的 ODE 初值问题；需要中高精度轨迹；Lotka–Volterra、SIR、机械系统、经济动态等光滑右端；需要事件检测（阈值穿越）、非负约束与密集输出。
核心优点	RK4 固定步长时全局误差 $O(h^4)$ ；`ode45` 自动步长，使用省心；单步法，易起步、易变步长；生态成熟，文档与例程丰富。
局限与缺点	显式 RK 对刚性效率差；每步多次计算 $f$ （RK4 为 4 次）；不是为守恒结构专门设计的辛方法；容差过松时仍可能在峰值处失真。

使用前应确认：

已把模型写成一阶系统 $\mathbf{y}'=\mathbf{f}(t,\mathbf{y})$ ；
粗判是否刚性（快慢时间尺度比是否极大）；
设定有意义的 RelTol/AbsTol，而不是完全放任默认而不检查；
准备守恒量、解析极限或容差加密验证。

三、核心原理：Butcher 表与自适应

1. 经典 RK4

\begin{aligned} k_1 &= f(t_n,y_n),\\ k_2 &= f\!\left(t_n+\tfrac{h}{2},\, y_n+\tfrac{h}{2}k_1\right),\\ k_3 &= f\!\left(t_n+\tfrac{h}{2},\, y_n+\tfrac{h}{2}k_2\right),\\ k_4 &= f\!\left(t_n+h,\, y_n+h k_3\right),\\ y_{n+1} &= y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4). \end{aligned}

局部截断误差 $O(h^5)$ ，全局误差 $O(h^4)$ 。

2. 一般 s 级 RK

\begin{aligned} k_i &= f\!\left(t_n+c_i h,\, y_n+h\sum_{j=1}^{s}a_{ij}k_j\right),\\ y_{n+1} &= y_n+h\sum_{i=1}^{s}b_i k_i. \end{aligned}

系数 $(a_{ij},b_i,c_i)$ 构成 Butcher 表。显式方法要求 $j\ge i$ 时 $a_{ij}=0$ ，从而可顺序计算。

3. 自适应步长思想

若同一步有 $p$ 阶与 $p+1$ 阶两个结果 $y$ 与 $\hat y$ ，则

\mathrm{err}\approx\|\hat y-y\|

可作为局部误差指示。若 $\mathrm{err}$ 大于容差则缩小 $h$ 重算，否则接受并可能放大下一步 $h$ 。ode45 即按此工作。

4. 与多步法对比

项目	RK（单步）	Adams（多步）
启动	只需当前点	需要多个历史点
变步长	自然	系数需重建
每步 $f$ 次数	较多	预报校正后可较少
默认工具	`ode45`	`ode113`

四、关键决策

求解器选择
- 非刚性：ode45（默认）
- 可能刚性：ode15s
- 高精度非刚性：ode113
- 必须展示公式：自写 RK4
容差
odeset('RelTol',1e-6,'AbsTol',1e-8) 是常用起点。容差收紧 10 倍后关键指标变化很小时，可认为数值误差可控。
选项技巧
- NonNegative：种群、仓室变量；
- Events：峰值时间、阈值穿越；
- MaxStep：强制不超过观测采样间隔；
- Mass：微分代数问题（进阶）。
输出
指定 tspan 为两点让求解器自适应内部步，再用 deval 或把 tspan 设为输出网格；不要误以为输出点就是全部内部步。
验证
至少做：容差加密、守恒量检查、与 RK4 固定步长或另一求解器交叉。

五、MATLAB 示例

%% 1) 自写 RK4：y'=-y
f = @(t, y) -y;
t0 = 0; tEnd = 2; y0 = 1; h = 0.1;
t = t0:h:tEnd;
y = zeros(size(t)); y(1) = y0;
for n = 1:numel(t)-1
    tn = t(n); yn = y(n);
    k1 = f(tn, yn);
    k2 = f(tn + h/2, yn + h/2*k1);
    k3 = f(tn + h/2, yn + h/2*k2);
    k4 = f(tn + h,   yn + h*k3);
    y(n+1) = yn + h/6*(k1 + 2*k2 + 2*k3 + k4);
end

figure;
plot(t, y, "o", t, exp(-t), "-"); grid on;
legend("RK4", "解析解");
title(sprintf("RK4 最大误差 = %.3e", max(abs(y - exp(-t)))));

%% 2) ode45 求解阻尼非线性摆
b = 0.25; g = 9.81; L = 1;
fPend = @(t, z) [ z(2);
    -(b)*z(2) - (g/L)*sin(z(1)) ];
opts = odeset("RelTol", 1e-7, "AbsTol", 1e-9);
[t, z] = ode45(fPend, [0 20], [pi/2; 0], opts);

figure;
subplot(2,1,1);
plot(t, z(:,1), "LineWidth", 1.2); grid on;
ylabel("\theta(t)"); title("ode45：阻尼摆角度");
subplot(2,1,2);
plot(z(:,1), z(:,2), "LineWidth", 1.2); grid on; axis equal;
xlabel("\theta"); ylabel("\omega"); title("相轨线");

%% 3) 容差收敛：SIR 峰值感染人数
N = 1e5; beta = 0.4; gamma = 0.1;
fSIR = @(t, y) [ -beta*y(1)*y(2)/N;
                  beta*y(1)*y(2)/N - gamma*y(2);
                  gamma*y(2) ];
y0 = [N-10; 10; 0];
relTols = [1e-3, 1e-4, 1e-5, 1e-6];
Ipeak = zeros(size(relTols));
for i = 1:numel(relTols)
    op = odeset("RelTol", relTols(i), "AbsTol", 1e-8, ...
        "NonNegative", 1:3);
    [~, Y] = ode45(fSIR, [0 200], y0, op);
    Ipeak(i) = max(Y(:,2));
end
disp(table(relTols(:), Ipeak(:), 'VariableNames', {'RelTol', 'I_peak'}));

%% 4) 事件定位：感染人数首次回落到 50
function [value, isterminal, direction] = eventI(t, y)
value = y(2) - 50;   % I - 50 = 0
isterminal = 0;      % 不终止
direction = -1;      % 从上往下穿越
end

opE = odeset("RelTol", 1e-6, "Events", @eventI, "NonNegative", 1:3);
[t, Y, te, Ye] = ode45(fSIR, [0 200], y0, opE);
disp("I 下行穿越 50 的时刻：");
disp(te');

RK4 封装：

function [t, Y] = rk4(f, tSpan, y0, h)
t = tSpan(1):h:tSpan(2);
if t(end) < tSpan(2) - 10*eps, t = [t, tSpan(2)]; end
Y = zeros(numel(t), numel(y0));
Y(1,:) = y0(:).';
for n = 1:numel(t)-1
    hn = t(n+1) - t(n);
    yn = Y(n,:).'; tn = t(n);
    k1 = f(tn, yn);
    k2 = f(tn+hn/2, yn+hn/2*k1);
    k3 = f(tn+hn/2, yn+hn/2*k2);
    k4 = f(tn+hn,   yn+hn*k3);
    Y(n+1,:) = (yn + hn/6*(k1+2*k2+2*k3+k4)).';
end
end

六、如何解释结果

使用 RK/ode45 时，规范表述包括：

求解器名称与阶数思想（RK4 或 Dormand–Prince 5(4)）；
RelTol、AbsTol、是否非负、是否事件；
积分区间、初值、关键输出指标；
容差加密表：证明数值误差小于建模结论所需精度；
与机制相关的解读（峰值、终态、振荡衰减率）；
若改用 ode15s 后结果变化大，讨论刚性可能。

注意：ode45 返回成功并不自动验证模型。应单独检查仓室和、能量、边界行为等。

七、常见误区

1. 默认容差一跑了之

默认值对许多题目够用，但对尖锐峰值或参数辨识可能不够。关键论文结论应附容差敏感性。

2. 刚性系统死用 `ode45`

步数爆炸或极慢时，优先怀疑刚性，换 ode15s 并比较结果。

3. 把输出网格当成内部步长

稀疏输出点不代表求解器用了大步长；内部步由误差控制决定。

4. 自写 RK4 时 $k_i$ 共用了被覆盖的 $y$

必须用临时变量按公式递推，避免原地修改导致阶数退化。

5. 认为高阶方法可以无视稳定性

显式 RK 仍有稳定域；大 Lipschitz 下过大固定 $h$ 仍会炸。

八、小结

Runge–Kutta 法的核心思想可以概括为：

在一个时间步内对右端函数进行多次取样并加权，用更高阶的方式逼近状态增量；自适应 RK 则进一步用不同阶公式的差来控制步长。

对数学建模中的绝大多数非刚性 ODE，ode45 是合理默认；自写 RK4 适合教学、对照与完全可控的实现。无论哪种，都要把容差、守恒与机制验证写进论证，而不是只展示一条光滑曲线。

一、直观理解：一步之内多点取样再加权​

二、适用场景与优缺点​

三、核心原理：Butcher 表与自适应​

1. 经典 RK4​

2. 一般 s 级 RK​

3. 自适应步长思想​

4. 与多步法对比​

四、关键决策​

五、MATLAB 示例​

六、如何解释结果​

七、常见误区​

1. 默认容差一跑了之​

2. 刚性系统死用 ode45​

3. 把输出网格当成内部步长​

4. 自写 RK4 时 kik_iki​ 共用了被覆盖的 yyy​

5. 认为高阶方法可以无视稳定性​

八、小结​