简介：自适应变步长的龙格库塔法是数值积分与常微分方程求解中的常用算法，这份MATLAB代码包将核心思路整理为可直接运行的脚本和说明，适合正在学习数值分析、需要将理论转换为程序实现的开发者参考。包体非常小巧，共4个文件，包括3个.m函数/脚本和1个txt程序说明，压缩包仅2KB，便于快速阅读和移植。代码覆盖了龙格库塔公式、步长动态调整与误差估计等关键环节，结合说明文件可理清从定义方程到主循环的整体结构。已有989人学习下载，对于想掌握自适应步长控制策略或开展相关课程设计的用户，是一份轻量且针对性强的入门范例。
1. 自适应变步长的龙格库塔法：为什么固定步长反而更难用
如果你调过ode45，其实已经用上了自适应变步长的龙格库塔法。MATLAB 里绝大多数求解器都不是按固定步长推进的，而是每一步都估算局部误差，误差太大就把h缩小重算，误差太小就把h放大省时间。标题里这个“自适应变步长的龙格库塔法”指的就是这类算法的最朴素实现：用定步长 RK4 做底子，自己写一套步长控制逻辑，最终得到一个和ode45行为类似、但每一步由你自己掌控的积分器。适合的场景很明确：你不想被ode45的封装接口限制，想在步长序列、误差估计、函数调用次数上拿到细粒度数据，或者在纯 MATLAB 课程设计里演示误差控制和收敛性分析。这篇文章先把“误差怎么估、步长怎么改”这两个核心问题讲透，再给一组可以直接复制运行的 MATLAB 函数，最后用ode45做基准对比和异常排查。
2. 从定步长 RK4 到自适应变步长：步长控制的误差估计原理
2.1 为什么定步长 RK4 在较长积分区间上会两头吃亏
经典四阶龙格库塔法做一步的公式是标准的四级四阶格式，单步精度高但没有任何误差提示。实际使用时只能提前猜测一个全局统一的h：取小了，积分区间稍微拉长，迭代次数就上千，数值舍入误差也会累积；取大了，在一些解变化剧烈的区段，比如 Van der Pol 方程的快速上升沿，RK4 会直接越过陡峭变化，生成明显错误的曲线，而且你很难从结果上看出问题。
自适应变步长的核心思路不是“每一步都精确”，而是把每一步的局部截断误差限制在用户给的容差范围内。每一步算出误差估计后，如果误差偏大就拒绝这一步并缩小步长，如果误差远小于容差就可以放大步长，用更少的步数完成整个区间。这样做的直接好处是：积分器把计算资源集中到解变化最剧烈的地方，而在平缓区段用大步长提高效率。一个工程上常见的反直觉结论是，自适应变步长代码的 CPU 时间往往少于定步长 RK4，因为平缓区段省下的步数远超陡峭区段增加的步数。
2.2 用步长折半的差分结构估计局部截断误差
要实现自适应，首先要回答一个问题：怎么知道这一步算得准不准？最常见的做法是步长折半差分。设当前状态为y_n，时间步长为h，先用完整步长h做一次 RK4，得到y_{n+1}^{(h)}；再把h拆成两个h/2，从y_n连续走两步，得到y_{n+1}^{(h/2)}。由于 RK4 的全局累积误差是O(h^4)，单步局部截断误差是O(h^5)，所以这两个结果之间的差值近似等于误差的（2^4 - 1）倍。
具体到代码里，局部误差估计写作：
lte = norm(y1 - y2, inf) / 15;这里y1是大步长结果，y2是两步半步长的结果，除以 15 是把二阶差分还原为一步误差。除以 15 这个细节很多人会漏掉，但它在步长控制里非常关键：如果不除，误差估计会被放大一个数量级，导致步长控制过度保守，最终积分步数激增。步长越小误差越小，折半两次得到的差值已经接近误差真值，这种估计方法在工程上足够用，误差量级正确，实现又最简单。
2.3 步长更新公式与安全因子的取值
有了局部误差估计err，下一步就是决定接受还是拒绝当前步。设用户给定的相对容差为RelTol，绝对容差为AbsTol，那么每一步的归一化误差可以写成：
sc = AbsTol + RelTol * max(abs(y), abs(y1)); err = max(abs(y1 - y2) ./ sc) / 15;当err <= 1时接受当前步，进入到下一步；否则拒绝，把h缩小重算。h的更新公式沿用控制器思想：
h_next = h * min(4, max(0.2, 0.9 / err^(1/5)));指数取1/5是因为 RK4 局部截断误差阶数是 5，要让err与h^5成正比。安全因子 0.9 是给误差估计的不确定性留余量，防止步长反复在拒绝和接受之间振荡；上限 4 和下限 0.2 则是防止单次步长突变过大。下面这条曲线行为在调试时值得记住：如果err经常恰好落在 0.9 到 1.1 附近，说明安全因子偏小或容差设置太紧，应该先调安全因子而不是调容差。
3. 在 MATLAB 里写出可运行的自适应变步长龙格库塔函数
3.1 函数签名与模块划分设计
先确定接口。这个函数要能控制初始步长，能返回每一步的时间点和状态，最好还能返回函数调用次数，方便对比效率。我通常会这样设计签名：
function [t, y, info] = adaptive_rk4(f, tspan, y0, RelTol, AbsTol, h0)其中f是微分方程句柄，格式为dy = f(t, y)；tspan是[t0, tf]；y0是初始状态列向量；RelTol和AbsTol是容差，h0是初始步长。info是结构体，里面记录fCalls（函数调用次数）和steps（接受步数）。内部实现拆成两层：主循环负责步长控制，局部函数rk4_step负责单步推进。这样主循环的代码量小，出问题时定位快。
MATLAB 里每一步要调用两次rk4_step，一次大步长、一次两个半步长。稍微算一下就知道，一次完整试算需要调用 4 次大步步进函数、8 次半步步进函数，也就是 12 个右端项函数值，而ode45的 Dormand-Prince 对只需要 6 个右端项就能同时得到结果和误差估计。这就是为什么工程上ode45用嵌入格式而不是步长折半法。这里自己写步长折半法，价值不在性能，而在于结构简单、每一步在做什么完全透明。
3.2 核心循环代码及参数说明
下面是完整的自适应变步长 RK4 实现，严格区分变量名，方便逐行核对误差估计过程：
function [t, y, info] = adaptive_rk4(f, tspan, y0, RelTol, AbsTol, h0) % 基于定步长RK4 + 步长折半误差估计的自适应积分器 t0 = tspan(1); tf = tspan(2); t = t0; y = y0(:); h = h0; % 控制参数 fac = 0.9; % 安全因子 facmin = 0.2; % 步长缩小下限倍数 facmax = 4.0; % 步长放大上限倍数 order = 4; fCalls = 0; while t < tf % 最后一步不能越过积分终点 h = min(h, tf - t); accepted = false; while ~accepted % 大步长：从 t 出发走 h [y1, ~, n1] = rk4_step(f, t, y, h); % 两个半大步长：先走 h/2，再走 h/2 [ym, ~, n2] = rk4_step(f, t, y, h/2); [y2, ~, n3] = rk4_step(f, t + h/2, ym, h/2); fCalls = fCalls + n1 + n2 + n3; % 归一化误差估计，除以 15 是关键 sc = AbsTol + RelTol * max(abs(y), abs(y1)); err = max(abs(y2 - y1) ./ sc) / 15; if err <= 1 % 接受当前步 t = t + h; y = y1; accepted = true; end % 更新步长：无论接受与否都按误差重新计算 if err > 0 h = h * min(facmax, max(facmin, fac / err^(1/(order+1)))); else h = h * facmax; % 误差为0，直接放大 end end end info.fCalls = fCalls; info.steps = numel(t) - 1; end function [yend, K, nf] = rk4_step(f, t, y, h) % 定步长四阶龙格库塔单步推进 k1 = f(t, y); k2 = f(t + h/2, y + h*k1/2); k3 = f(t + h/2, y + h*k2/2); k4 = f(t + h, y + h*k3); yend = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6; K = [k1, k2, k3, k4]; nf = 4; end这段代码的逻辑说明：主循环先做一个完整大步长、两个半大步长，用三次rk4_step调用得到两个候选解；归一化误差err是一个标量，由所有状态分量中最大的相对偏差决定，这是norm(..., inf)的语义，它比均方根更严格。步长更新在if之外，所以即使当前步被拒绝，新的h也已经算出，不需要再做一次求幂运算。注意y1和y2如果某个分量恰好穿过零点，max(abs(y), abs(y1))可以避免sc太小导致误差被放大；绝对容差AbsTol默认给 1e-6 量级，如果被积状态本身很小，要相应调小。
3.3 对输出节点的重采样处理
上面的函数返回的是所有被接受步的节点，但这些节点分布不均匀，做图时如果直接连线，曲线在陡峭区段会显得比其他区段“密”，这符合物理规律，但如果你要求等间隔输出，比如控制周期固定为 0.01 秒，就必须做后处理。最简单的方法是在主循环尾部加一个输出插值需求：先记录t_raw、y_raw，积分结束后用interp1(t_raw, y_raw, t_query)在等间隔时间点取值。要注意的是interp1要求节点单调递增，自适应步长数组天然满足这个条件；另外插值精度和一阶线性插值一致，若需要高阶精度，可以把spline打开，但这会引入轻微的过冲，在斜坡类响应上要小心。
不想自写插值也可以换一种做法：直接把tspan拆成两段，比如先积到[t0, t_mid]，再积到[t_mid, tf]，最终把两组输出拼接。这种做法每段的终值来自自适应积分精度保证，拼接处不会像插值那样产生额外误差，但会让平缓区段被迫加密，失去自适应省步数的意义。常见工程做法是保留adaptive_rk4的原始输出，仅在需要展示标准间隔曲线时调用一次interp1。
4. 用 ode45 对比验证自适应变步长龙格库塔代码，并排查 error 9 类异常
4.1 与 ode45 的精度和调用次数对照
写完代码后，第一步不是直接换到自己的项目里，而是拿一个已知解析解的非线性方程做对照。我拿下面的 Logistic 增长方程做测试：
dy/dt = 0.1 * y * (1 - y/100), y(0) = 1, 0 <= t <= 200这个方程解析解是 S 型曲线，积分区间跨越快速增长段和饱和段，非常适合检验自适应步长是否真的在陡峭区段加密。分别调用adaptive_rk4和ode45，两者都设RelTol = 1e-6, AbsTol = 1e-8，统计末端值误差和函数调用次数。典型的对照结果落在下面这个区间内：
| 求解方式 | 函数调用次数占比 | 末端相对误差 | 步长变化范围 |
|---|---|---|---|
| ode45 (DOPRI5) | 基准 | 约 1e-7 量级 | 自动 |
| 步长折半自适应 RK4 | 约为 ode45 的 1.5 到 2 倍 | 约 1e-7 到 1e-6 量级 | 约为 ode45 的 0.5 到 2 倍 |
| 定步长 RK4, h=0.1 | 接近自适应 RK4 | 可能超过 1e-4 | 固定 |
函数调用次数多出一半左右是符合预期的，因为步长折半需要 12 次右端项计算，而ode45嵌入对只需要 6 次。注意末端相对误差仍然在容差范围内，说明步长折半的误差估计机制有效。这个对比的价值在于：当你想替换ode45时，先知道自己的代码要多花多少函数调用，并设置一个可接受的预算，后面做实时仿真才有的放矢。
4.2 状态不连续导致 NaN 传播的典型路径
自适应步长代码在工程里最常见的异常不是“报错”，而是静默生成NaN或Inf。首步正常，从某一步开始y全部变成NaN，然后步长控制因为err也是NaN而走进死循环。造成这种问题的最常见原因是右端项函数存在不连续，比如碰撞、死区、离散事件，而自适应步长在事件前后仍然按连续系统控制误差。误差估计会突然暴涨，步长被迫缩到极小，极端情况下err变成Inf，1 / err^(1/5)直接变 0，步长被压到浮点数下溢。
排查顺序给到这里：先检查f里有没有除零、有没有对负数开偶次方、有没有log(0)；再把adaptive_rk4内部每一步的h、err打印出来，定位从哪一步开始出现非有限值；最后将事件附近的时间点打印出来，确认不连续性是否发生在t的某个精确值附近。如果确实有事件，用ode45的Events选项把积分在事件点停下来，处理完事件后重新启动积分器，这比在右端项函数里硬编码if语句更可控。
至于终端偶发出现的 Error 9 类报错，常见于 MEX 文件或文件 ID 操作场景，如果你没改info里的文件句柄，多半是系统层面的文件描述符问题。先确认 MATLAB 工作目录和项目路径是否可写，再查是否打开了大量文件没关闭。这类报错通常和积分算法无关，却被误认为是自适应步长的 bug，浪费不少调试时间。
4.3 刚性问题和 ode15s 的选用边界
自适应步长 RK4 在刚性方程上会彻底失效。判断方法很直接：记录adaptive_rk4输出的最小步长h_min和积分区间长度tf - t0，如果h_min在1e-8量级而tf是10，说明求解器在把步长压到荒谬的小。像 Van der Pol 在mu = 1000的振荡、化学动力学里的快慢反应耦合，都属于典型的刚性系统。这类问题要用隐式方法，常见选择是ode15s或ode23t。
当你的项目从自适应变步长 RK4 切换到ode15s时，参数习惯要改：RelTol默认保持1e-3即可，没必要追求1e-6，因为隐式方法每一步要解线性方程组，RelTol每缩小一个数量级，Jacobian 计算和分解的成本会明显上涨。自适应步长 RK4 适合的是非刚性、中等精度、想看清步长行为的教学代码和轻量仿真。如果你发现自己不断降低容差来压制振荡，先检查模型是不是刚性，而不是继续调 RK4。
5. 自适应变步长龙格库塔代码的调优技巧：批量评估与容差手感
5.1 用批量状态评估代替单点循环
写自适应积分器时，右端项函数是最容易被拖慢的瓶颈。很多人直接把f写成只接受单列向量的函数，主循环里一步调 12 次，一个 5000 步的积分会触发几万次函数调用。更快的常见做法是在rk4_step里允许y按列拼接成矩阵，一次调用同时评估多个状态候选值。修改f的签名就可以做到批量评估，例如F = f(t, [y, y + h*k1/2, y + h*k2/2, y + h*k3])，一次右端项调用得到四个列向量，再把它们按 RK4 系数线性组合。这种改法在 MATLAB 里能利用内置向量化加速，比for循环包四列快不少，代价是代码可读性下降。建议保留一个单点评估版本用于调试，批量版本用于长区间积分。
5.2 容差设置的经验区间
容差不是越小越好。RelTol从 1e-3 调到 1e-6，误差会下降，但步长变小，函数调用次数上升。实际经验是：画趋势图用RelTol=1e-3，做定量分析用1e-6，高于1e-9的情况极少，除非你在做高精度轨道外推。绝对容差AbsTol要匹配状态量纲，如果状态里有幅值为 1e-5 的分量，AbsTol=1e-8会迫使积分器做很多无效小步长，这种情况应该把AbsTol降到1e-10，甚至给状态向量里每个分量单独设置容差。我的习惯是把AbsTol写成向量，与y0对齐，这样不同量纲的状态各得其所。
5.3 最后做一次“步长曲线目检”
调完容差和批量评估后，把adaptive_rk4输出的步长序列画出来。横轴是时间，纵轴是步长h。一条健康的步长曲线应该是：平缓区段步长大，陡峭区段步长小，两者过渡平滑，没有频繁的锯齿。如果步长在每个积分点上都上下跳动，说明安全因子 0.9 太小或误差估计有偏，先调大安全因子到 0.95 看看是否稳定；如果步长在某一段持续保持在facmin下限，说明该区段可能存在不连续点或模型刚性。这个目检操作 30 秒就能完成，却是判断自适应步长龙格库塔代码质量最直接的验证手段，比盯着末端误差一个数字有效得多。
本文还有配套的精品资源，点击获取
