简介一份面向MATLAB数值计算学习者的自适应变步长龙格库塔法源码包用于解决常微分方程初值问题中固定步长算法精度与效率难以兼顾的问题。资源核心是完整的MATLAB实现包含方程定义、龙格库塔公式、步长自适应调节和主程序调用等模块适合数值分析课程实践或科研仿真中快速嵌入使用。压缩包内共4个文件包含3个MATLAB源文件.m和1个txt程序说明源文件涵盖算法主程序、函数定义与步长控制逻辑说明文件便于快速了解调用方式和参数设置整包仅2KB轻量精炼适合按需修改。已有989人学习下载代码同时演示了误差估计和步长调整策略读者既可以获得能直接运行的算法模板也能掌握自适应步长控制的实现思路并迁移到电路分析、流体模拟等实际仿真场景。1. 自适应变步长的龙格库塔法当同一套 RK4 需要对付一百倍速度差异自适应变步长的龙格库塔法看起来只是在固定步长前面加了一个“自适应”前缀实际跑起来差别极大。用 MATLAB 做数值仿真时很多方程同时包含快速瞬态和平缓段比如 t0 附近需要 h1e-4 才能抓住曲线转折进入平缓段后 h 取 0.1 都绰绰有余。若固定步长要么前期无法收敛要么后期做大量无效计算。这个压缩包里的 zhukong.m 是主控循环half.m 提供半步长推进f.m 只放右端项程序说明.txt 则把这套配合关系做了文字交代。对于刚接触自适应积分、想看清步长控制逻辑、又不想直接调 ode45 黑盒的从业者来说这套代码是一份能下断点的白盒样本。2. 从 RK4 的局部误差看步长为什么能自动调整2.1 龙格库塔法的递推骨架与 Butcher 表自适应变步长的基本判断依据来自同一套 RK4 递推骨架。固定步长的 RK4 更新公式写成 MATLAB 语法的形式k1 f(t, y); k2 f(t h/2, y h/2*k1); k3 f(t h/2, y h/2*k2); k4 f(t h, y h*k3); y_next y h/6*(k1 2*k2 2*k3 k4);这里的 k1 到 k4 并不是随意取四个点做平均而是把区间 [t, th] 内部几个典型位置的切线方向按不同权重组合起来。一阶欧拉法只使用起点斜率大步长时方向误差直接累积RK4 通过四次求值把泰勒展开匹配到四阶局部截断误差项是 O(h^5)。这个五阶项就是自适应步长做判断的“原料”。细化到系数层面可以用 Butcher 表表达同一组系数0 | 1/2 | 1/2 1/2 | 0 1/2 1 | 0 0 1 ---------------------- | 1/6 1/3 1/3 1/6表左列是 c_i表示每个 k_i 在时间轴上的位置中间矩阵 a_ij 表示上一个 k 对当前 k 的线性组合底部一行 b_i 是输出加权系数。定步长程序只需要按这张表从上到下算出 k_i最后做加权和。自适应程序的不同在于它要多算一条参考线用两条不同路径的结果相减得到误差估计。2.2 全步长减半步长误差怎么算出来误差估计不一定要用更高阶公式同一个 RK4 也能通过“走的路径不同”制造出差值。常见做法是在 [t_n, t_nh] 上直接算一步完整 RK4得到 y_full再把区间切成两段每段 h/2 各做一次 RK4得到 y_half。以这两个结果之差 err abs(y_full - y_half) 作为局部误差信号。为什么这个差值能代表误差记真实解为 y_exactRK4 的局部误差是 O(h^5)所以y_full y_exact C * h^5 O(h^6) y_half y_exact C * (h/2)^5 O(h^6)这里 C 是某个和方程导数有关的常数不必要知道具体值。两式相减后主项是 (31/32) C h^5它和真正的误差 C h^5 同量级。程序把这个差值当作“误差超标”的指标差值小说明两条路径给出的解接近步长可以放大差值大说明当前区间内解的变化剧烈步长必须缩小。需要提醒一个容易踩的地方误差范数的选择。上面写成 abs 只适合标量 y对于方程组应该使用带权范数而不是对每个分量单独判断。否则会出现某个分量误差很大但因为它和其他分量量级不同abs 判断完全不敏感的情况。一般做法是err sqrt(sum(((y_full - y_half) ./ (atol rtol*abs(y_half))).^2));atol 是绝对容限rtol 是相对容限。这样误差信号的含义是“相对解尺度的误差”而不是纯粹数值差。对量级跨越很多个数量级的模型只用 abs(y_full - y_half) 会导致小量级段始终判断误差超限步长被反复压缩计算效率急转直下。2.3 步长调整公式与安全系数有了 err 以后把 err 和容限阈值 1归一化后比较err 1 接受当前步并放大 herr 1 拒绝当前步并缩小 h。缩放指数从误差阶来RK4 的局部误差与 h^5 成正比所以让步长乘上一个系数 r误差大约变成原来的 r^5。要满足 r^5 * err 1则 r err^(-1/5)。实际代码为了避免过激调整会加安全系数h_new h * min(max(0.9 * err^(-1/5), 0.2), 5.0);这个式子里的 0.9 是经验安全系数下标 0.2 和上标 5.0 把单步调整限制在一个可控范围。下面这张表总结了最常见的参数选择参数典型值作用tol1e-6单步误差阈值越小越逼近精确解err 计算方式abs 或加权范数决定误差信号代表绝对差还是相对差指数 1/5与 RK4 局部误差阶数对应指数写错会导致步长过冲或缩过头放大上限5.0防止平缓段步长暴涨越过快速变化区缩小下限0.2避免单步把步长缩得太狠等待多步连续收缩安全系数0.8~0.9给误差留裕量减少接受/拒绝的来回振荡这个指数 1/5 是最容易被改成 1/4 的地方。改小后放大步长时步子迈得不够收缩步长时又缩不过头可能让整段积分步数多出几倍。后面排错时如果发现 h 序列出现抖动先检查是不是这里写成了 1/4。3. zhukong.m、half.m、f.m三件套如何拼出一个自适应循环3.1 f.m主循环只认识“求导”这一件事在压缩包的三文件结构里f.m 是唯一和具体方程打交道的地方。它不用关心步长也不用关心误差只需要返回右端项的导数值function dydt f(t, y) % 示例方程刚性衰减到 cos(t) dydt -1000 * (y - cos(t)) - sin(t); end说明t 是当前时刻y 是状态变量。这里写的是一个标量方程但建议把 y 设计成列向量因为 half.m 和 zhukong.m 的主循环会把 y 当作状态向量做线性组合。方程右侧带 -1000 的系数是用来测试自适应的典型刚性问题t0 附近解快速跌落到 cos(t)之后基本跟着 cos(t) 平缓波动。这样一条方程可以直观看到步长从 1e-4 级别涨到 0.1 级别的全过程非常适合当白盒测试用例。f.m 不要内置任何全局变量或打印语句。自适应循环可能在同一时刻反复调用 f刚进入时计算新步长失败回退后又要重算如果 f 里有输出语句控制台会被刷爆。程序说明.txt 如果要求改模型只需要改这一个文件的返回值。3.2 half.m半步长 RK4 推进函数half.m 名字来自“半步长”这个策略。它把 RK4 单步封装成一个函数专门负责从任意起点 t0 推进到 t0h内部只有四次 f 调用function y_next half(t0, y0, h, f) % 对 [t0, t0h] 做一次 RK4 推进 k1 f(t0, y0); k2 f(t0 h/2, y0 h/2*k1); k3 f(t0 h/2, y0 h/2*k2); k4 f(t0 h, y0 h*k3); y_next y0 h/6*(k1 2*k2 2*k3 k4); end这个函数只做一件事从 y0 推进到 t0h。它不知道外面的主循环在计算全步长还是半步长也不知道误差是否超限。参数 t0 是当前时刻y0 是当前状态列向量h 是这段区间的长度f 是函数句柄。类型上要保证 f(t, y) 返回和 y 同形状的向量否则 y_next 的行列维度会出问题。在 zhukong.m 里这个函数会被调用三次一次算全步长 h一次算前半段 h/2一次算后半段 h/2。之所以不把第三个调用合并是为了让全步长路径和半步长路径在算法上尽量独立避免复用中间 k 值改变误差估计的统计特性。虽然复用 k1 可以省一次 f 调用但它会让两条路径的误差信号失去独立性接近阈值时更容易产生误判。3.3 zhukong.m接受-拒绝循环是真正的控制器zhukong.m 是整套代码的控制中心它调用 half.m 获得 y_full 和 y_half再根据误差信号决定接受还是回退。一个可以放进 MATLAB 运行的最小主循环如下function [T, Y, H] zhukong(f, tspan, y0, h0, tol) T(1) tspan(1); Y(1,:) y0(:); H(1) h0; t tspan(1); y y0(:); h h0; maxit 200000; for k 1:maxit if t tspan(2) break; end if t h tspan(2) h tspan(2) - t; end y_full half(t, y, h, f); y_mid half(t, y, h/2, f); y_half half(t h/2, y_mid, h/2, f); err abs(y_full - y_half) / max(tol, tol*abs(y_half)); if err 1 t t h; y y_half; T(end1,:) t; Y(end1,:) y(:); H(end1,:) h; h h * min(max(0.9 * err^(-1/5), 0.2), 5.0); if err 1e-8 % 误差远小于阈值说明还有放大余量 h h * 1.5; end else h h * max(0.2, 0.9 * (1 / err)^(1/5)); end end end这段代码的逻辑可以用接受/拒绝两条路径归纳。接受路径把 y_half 作为当前解并放大步长拒绝路径不推进 t只缩小 h然后重新计算逻辑回到同一个 t。需要注意 err 这行是标量写法针对方程组需要换成加权范数 norm(...)否则 MATLAB 不会进入 if 分支。变量含义和需要警惕的点列成一张表变量含义需要警惕的点h0初始步长不要取得比物理时间尺度大很多否则第一轮就拒绝tol归一化容限tol 太大会把明显误差当可接受y_full / y_half全步长和半步长结果接受时必须用 y_half而不是更粗糙的 y_fullmaxit循环上限刚性方程拒绝次数多maxit 要留足H 输出每步实际使用的步长用于事后画步长序列定位刚区间额外的 1.5 倍放大分支不要滥用。上面代码在 err 1e-8 时又乘了 1.5这个分支在病态问题里可能让 h 从刚缩小的状态突然反弹导致下一次误差重新超限。实际使用时我更倾向去掉这个分支只保留 err^(-1/5) 的 0.9 倍系数让步长自然增长。4. 容限、初始步长和缩放因子边界设错了照样跑不动4.1 绝对容限和相对容限怎么配合自适应变步长误差阈值 tol 在很多教材里只有一个数字但落到 MATLAB 里要区分绝对容限和相对容限。绝对容限 atol 决定解接近 0 时允许的绝对偏差相对容限 rtol 决定解量级较大时允许的相对偏差。例如在电路仿真里电压和电流两个状态量量级差上百倍如果只用单一 abs 误差电流分量会被误判成误差很大程序会不断把步长缩小最后卡死在初始化阶段。常见的处理是沿用 ode45 的参数习惯rtol 1e-6; atol 1e-8; err sqrt(sum(((y_half - y_full) ./ (atol rtol*abs(y_half))).^2));注意 atol 和 rtol 都是一个标量作用于所有分量这其实是一种简化。更细的用法是让 atol 和 rtol 与 y 同维度对不同分量分别设置容限。比如初值问题中某个状态只变化 1e-5 量级那就单独给它一个 atol1e-7而另一个状态变化到 1e3 量级仍用 rtol 控制。分量级容限能显著减少刚性方程组里的无效拒绝。4.2 初始步长 h0 的经验做法自适应步长理论上是自动找步长但初始步长给得离谱会让循环在开始阶段空转几轮。最稳妥的经验是先跑一步欧拉试探取一个很小的临时步长 h_try算出 dydt 的初值用 dydt 变化率反推初始 h。对罗森布罗克这类刚性求解器更常见的策略是 h0 sqrt(eps) * max(norm(y0), 1e-6)。对本项目这种基于 RK4 的代码我一般直接给h0 1e-3 * max(norm(y0), 1);这个值对大多数非病态问题能在前几步内快速放大到合理区间。如果第一步就反复拒绝则说明 h0 比可接受范围大太多此时不要手动把 h0 调小到离谱的程度而是先检查 f.m 里的时间尺度。比如方程系数是 1e6而 h0 仍从 1e-3 出发第一步 f(t,y) 就会让 k2 和 k4 相差好几个数量级误差信号失真后面调整等于在噪声上做控制。更合理的做法是按 f 中最大时间常数的量级给初始步长。4.3 缩放系数、最小步长和最大步长的钳制步长更新公式 min(max(0.9 * (tol/err)^(1/5), low), high) 里的 low 和 high 不只是保护误差也保护程序本身。low 设置太小例如 1e-10一旦误差信号有异常峰值h 直接缩到接近机器精度后面就算误差恢复正常也要几百步才爬回来。更常见的做法是用 H 序列里最小步长做一个全局钳制当 h 1e-8 * tspan 时直接触发“求解失败”退出而不是继续死循环。if h 1e-8 * abs(tspan(2) - tspan(1)) warning(步长低于时间尺度的1e-8可能已经遇到奇点); break; end这里 1e-8 不是需要反复调节的神秘系数而是基于双精度浮点余量的保守选择。如果方程本身是刚性最小步长会在合法区间内自动降下来如果 h 降到这个阈值以下说明右端函数在这点附近有奇点或非连续再自适应也没有意义。high 上限方面不赞成把 high 设成 10 以上。自适应步长通过误差估计放大步长但如果放大速度过快区间内快速变化会被跨过去。比如 h 从 1e-4 放大到 1e-3放大 10 倍后误差也许刚好在阈值内但下一步误差突然超限程序又要缩回来形成一种低频振荡。上限 5.0 配合 0.9 的安全系数让 h 的增长路径更平滑。下面这张表可以作为参数调整入口调整目标改哪个参数副作用注意让结果更精确减小 atol / rtol总步数线性上涨刚性区间可能涨几倍让平缓段更快增大 high可能在快慢过渡区错过一个尖峰减少拒绝次数增大安全系数到 0.95单步误差变大但总调用次数可能下降解决启动失败改 h0 数量级只影响前几步不会改变后续步长轨迹最后要提醒在 zhukong.m 里错误地将最大步长 high 设成 100 并不会让程序跑得更快它只会让 t 推进的跳变变大然后每次误差检查都会失败最终总计算量反而增加。自适应方法的核心价值在于“该小的地方小该大的地方大”maxstep 上限只是一个防呆边界。5. 验算自适应行为用 H 输出序列反推刚区间位置5.1 把调参结果可视化zhukong.m 返回的 H 数组就是每步实际使用的步长。运行完求解后画两条线一条是解 y(t)一条是 h(t)刚区间立刻显形[T, Y, H] zhukong(f, [0 5], 1, 2e-4, 1e-6); subplot(2,1,1); plot(T, Y); title(y(t)); grid on; subplot(2,1,2); semilogy(T(2:end), H(2:end)); title(step size h(t)); grid on;这段代码用 semilogy 而不是 plot因为步长在最初几十步可能横跨 1e-5 到 1e-1对数纵轴才能看出指数级变化。如果发现 h(t) 在某个区间出现锯齿形骤降骤升说明误差在阈值附近来回震荡优先检查 err 那行的分母有没有用 atol rtol*abs(y_half)。5.2 与 ode45 做同点差值验证自适应实现写完后最直接的正确性验证是插值后和 ode45 结果对照。ode45 内部使用带误差控制的 Dormand-Prince 方法这套“全步长减半步长”是一对更朴素的策略两者不应该完全相同但应具有相近的数量级[T2, Y2] ode45(f, [0 5], 1, odeset(RelTol, 1e-6, AbsTol, 1e-8)); y_interp interp1(T, Y, T2, pchip); max(abs(y_interp - Y2))max 差值如果在 1e-4 量级以下说明主循环没有原则性错误。如果差值达到 1e-2第一步不是调小容限而是看 H 序列有没有大跳跃尤其关注 t0 附近前 10 步。多数实现问题都出在拒绝分支里忘记把 h 更新回下次迭代使用的变量导致实际步长和 H 记录不一致。5.3 把“拒绝次数”变成调参指标再给一个能直接量化的技巧在 zhukong.m 的 else 分支里加一个计数器统计整个求解过程中步长被拒绝的次数。这个数字比最终步数更能说明问题。容限设置合理时拒绝次数应远小于接受次数如果拒绝次数超过接受次数的 20%说明初始步长偏大或安全系数过小应当把 0.9 提到 0.95而不是把 tol 放大。相反如果 H 序列几乎不变说明步长一直被上限钳住此时应当调大 high或确认方程是否真的存在刚性区间。这一行放在 else 分支末尾即可reject_count reject_count 1;在函数返回前把 reject_count 作为第三个输出或写入全局变量。用静态分析看一个积分算法是好是坏单看总步数并不准确用“接受步数 / 总评估次数”来衡量更接近真实成本。half.m 每次调用会执行四次 f 求值一次全步长加两次半步长总计 12 次求值这个数才是和 ode45 对比性能时应关注的核心数值。本文还有配套的精品资源点击获取
自适应变步长龙格库塔法：从误差估计到MATLAB工程实现
简介一份面向MATLAB数值计算学习者的自适应变步长龙格库塔法源码包用于解决常微分方程初值问题中固定步长算法精度与效率难以兼顾的问题。资源核心是完整的MATLAB实现包含方程定义、龙格库塔公式、步长自适应调节和主程序调用等模块适合数值分析课程实践或科研仿真中快速嵌入使用。压缩包内共4个文件包含3个MATLAB源文件.m和1个txt程序说明源文件涵盖算法主程序、函数定义与步长控制逻辑说明文件便于快速了解调用方式和参数设置整包仅2KB轻量精炼适合按需修改。已有989人学习下载代码同时演示了误差估计和步长调整策略读者既可以获得能直接运行的算法模板也能掌握自适应步长控制的实现思路并迁移到电路分析、流体模拟等实际仿真场景。1. 自适应变步长的龙格库塔法当同一套 RK4 需要对付一百倍速度差异自适应变步长的龙格库塔法看起来只是在固定步长前面加了一个“自适应”前缀实际跑起来差别极大。用 MATLAB 做数值仿真时很多方程同时包含快速瞬态和平缓段比如 t0 附近需要 h1e-4 才能抓住曲线转折进入平缓段后 h 取 0.1 都绰绰有余。若固定步长要么前期无法收敛要么后期做大量无效计算。这个压缩包里的 zhukong.m 是主控循环half.m 提供半步长推进f.m 只放右端项程序说明.txt 则把这套配合关系做了文字交代。对于刚接触自适应积分、想看清步长控制逻辑、又不想直接调 ode45 黑盒的从业者来说这套代码是一份能下断点的白盒样本。2. 从 RK4 的局部误差看步长为什么能自动调整2.1 龙格库塔法的递推骨架与 Butcher 表自适应变步长的基本判断依据来自同一套 RK4 递推骨架。固定步长的 RK4 更新公式写成 MATLAB 语法的形式k1 f(t, y); k2 f(t h/2, y h/2*k1); k3 f(t h/2, y h/2*k2); k4 f(t h, y h*k3); y_next y h/6*(k1 2*k2 2*k3 k4);这里的 k1 到 k4 并不是随意取四个点做平均而是把区间 [t, th] 内部几个典型位置的切线方向按不同权重组合起来。一阶欧拉法只使用起点斜率大步长时方向误差直接累积RK4 通过四次求值把泰勒展开匹配到四阶局部截断误差项是 O(h^5)。这个五阶项就是自适应步长做判断的“原料”。细化到系数层面可以用 Butcher 表表达同一组系数0 | 1/2 | 1/2 1/2 | 0 1/2 1 | 0 0 1 ---------------------- | 1/6 1/3 1/3 1/6表左列是 c_i表示每个 k_i 在时间轴上的位置中间矩阵 a_ij 表示上一个 k 对当前 k 的线性组合底部一行 b_i 是输出加权系数。定步长程序只需要按这张表从上到下算出 k_i最后做加权和。自适应程序的不同在于它要多算一条参考线用两条不同路径的结果相减得到误差估计。2.2 全步长减半步长误差怎么算出来误差估计不一定要用更高阶公式同一个 RK4 也能通过“走的路径不同”制造出差值。常见做法是在 [t_n, t_nh] 上直接算一步完整 RK4得到 y_full再把区间切成两段每段 h/2 各做一次 RK4得到 y_half。以这两个结果之差 err abs(y_full - y_half) 作为局部误差信号。为什么这个差值能代表误差记真实解为 y_exactRK4 的局部误差是 O(h^5)所以y_full y_exact C * h^5 O(h^6) y_half y_exact C * (h/2)^5 O(h^6)这里 C 是某个和方程导数有关的常数不必要知道具体值。两式相减后主项是 (31/32) C h^5它和真正的误差 C h^5 同量级。程序把这个差值当作“误差超标”的指标差值小说明两条路径给出的解接近步长可以放大差值大说明当前区间内解的变化剧烈步长必须缩小。需要提醒一个容易踩的地方误差范数的选择。上面写成 abs 只适合标量 y对于方程组应该使用带权范数而不是对每个分量单独判断。否则会出现某个分量误差很大但因为它和其他分量量级不同abs 判断完全不敏感的情况。一般做法是err sqrt(sum(((y_full - y_half) ./ (atol rtol*abs(y_half))).^2));atol 是绝对容限rtol 是相对容限。这样误差信号的含义是“相对解尺度的误差”而不是纯粹数值差。对量级跨越很多个数量级的模型只用 abs(y_full - y_half) 会导致小量级段始终判断误差超限步长被反复压缩计算效率急转直下。2.3 步长调整公式与安全系数有了 err 以后把 err 和容限阈值 1归一化后比较err 1 接受当前步并放大 herr 1 拒绝当前步并缩小 h。缩放指数从误差阶来RK4 的局部误差与 h^5 成正比所以让步长乘上一个系数 r误差大约变成原来的 r^5。要满足 r^5 * err 1则 r err^(-1/5)。实际代码为了避免过激调整会加安全系数h_new h * min(max(0.9 * err^(-1/5), 0.2), 5.0);这个式子里的 0.9 是经验安全系数下标 0.2 和上标 5.0 把单步调整限制在一个可控范围。下面这张表总结了最常见的参数选择参数典型值作用tol1e-6单步误差阈值越小越逼近精确解err 计算方式abs 或加权范数决定误差信号代表绝对差还是相对差指数 1/5与 RK4 局部误差阶数对应指数写错会导致步长过冲或缩过头放大上限5.0防止平缓段步长暴涨越过快速变化区缩小下限0.2避免单步把步长缩得太狠等待多步连续收缩安全系数0.8~0.9给误差留裕量减少接受/拒绝的来回振荡这个指数 1/5 是最容易被改成 1/4 的地方。改小后放大步长时步子迈得不够收缩步长时又缩不过头可能让整段积分步数多出几倍。后面排错时如果发现 h 序列出现抖动先检查是不是这里写成了 1/4。3. zhukong.m、half.m、f.m三件套如何拼出一个自适应循环3.1 f.m主循环只认识“求导”这一件事在压缩包的三文件结构里f.m 是唯一和具体方程打交道的地方。它不用关心步长也不用关心误差只需要返回右端项的导数值function dydt f(t, y) % 示例方程刚性衰减到 cos(t) dydt -1000 * (y - cos(t)) - sin(t); end说明t 是当前时刻y 是状态变量。这里写的是一个标量方程但建议把 y 设计成列向量因为 half.m 和 zhukong.m 的主循环会把 y 当作状态向量做线性组合。方程右侧带 -1000 的系数是用来测试自适应的典型刚性问题t0 附近解快速跌落到 cos(t)之后基本跟着 cos(t) 平缓波动。这样一条方程可以直观看到步长从 1e-4 级别涨到 0.1 级别的全过程非常适合当白盒测试用例。f.m 不要内置任何全局变量或打印语句。自适应循环可能在同一时刻反复调用 f刚进入时计算新步长失败回退后又要重算如果 f 里有输出语句控制台会被刷爆。程序说明.txt 如果要求改模型只需要改这一个文件的返回值。3.2 half.m半步长 RK4 推进函数half.m 名字来自“半步长”这个策略。它把 RK4 单步封装成一个函数专门负责从任意起点 t0 推进到 t0h内部只有四次 f 调用function y_next half(t0, y0, h, f) % 对 [t0, t0h] 做一次 RK4 推进 k1 f(t0, y0); k2 f(t0 h/2, y0 h/2*k1); k3 f(t0 h/2, y0 h/2*k2); k4 f(t0 h, y0 h*k3); y_next y0 h/6*(k1 2*k2 2*k3 k4); end这个函数只做一件事从 y0 推进到 t0h。它不知道外面的主循环在计算全步长还是半步长也不知道误差是否超限。参数 t0 是当前时刻y0 是当前状态列向量h 是这段区间的长度f 是函数句柄。类型上要保证 f(t, y) 返回和 y 同形状的向量否则 y_next 的行列维度会出问题。在 zhukong.m 里这个函数会被调用三次一次算全步长 h一次算前半段 h/2一次算后半段 h/2。之所以不把第三个调用合并是为了让全步长路径和半步长路径在算法上尽量独立避免复用中间 k 值改变误差估计的统计特性。虽然复用 k1 可以省一次 f 调用但它会让两条路径的误差信号失去独立性接近阈值时更容易产生误判。3.3 zhukong.m接受-拒绝循环是真正的控制器zhukong.m 是整套代码的控制中心它调用 half.m 获得 y_full 和 y_half再根据误差信号决定接受还是回退。一个可以放进 MATLAB 运行的最小主循环如下function [T, Y, H] zhukong(f, tspan, y0, h0, tol) T(1) tspan(1); Y(1,:) y0(:); H(1) h0; t tspan(1); y y0(:); h h0; maxit 200000; for k 1:maxit if t tspan(2) break; end if t h tspan(2) h tspan(2) - t; end y_full half(t, y, h, f); y_mid half(t, y, h/2, f); y_half half(t h/2, y_mid, h/2, f); err abs(y_full - y_half) / max(tol, tol*abs(y_half)); if err 1 t t h; y y_half; T(end1,:) t; Y(end1,:) y(:); H(end1,:) h; h h * min(max(0.9 * err^(-1/5), 0.2), 5.0); if err 1e-8 % 误差远小于阈值说明还有放大余量 h h * 1.5; end else h h * max(0.2, 0.9 * (1 / err)^(1/5)); end end end这段代码的逻辑可以用接受/拒绝两条路径归纳。接受路径把 y_half 作为当前解并放大步长拒绝路径不推进 t只缩小 h然后重新计算逻辑回到同一个 t。需要注意 err 这行是标量写法针对方程组需要换成加权范数 norm(...)否则 MATLAB 不会进入 if 分支。变量含义和需要警惕的点列成一张表变量含义需要警惕的点h0初始步长不要取得比物理时间尺度大很多否则第一轮就拒绝tol归一化容限tol 太大会把明显误差当可接受y_full / y_half全步长和半步长结果接受时必须用 y_half而不是更粗糙的 y_fullmaxit循环上限刚性方程拒绝次数多maxit 要留足H 输出每步实际使用的步长用于事后画步长序列定位刚区间额外的 1.5 倍放大分支不要滥用。上面代码在 err 1e-8 时又乘了 1.5这个分支在病态问题里可能让 h 从刚缩小的状态突然反弹导致下一次误差重新超限。实际使用时我更倾向去掉这个分支只保留 err^(-1/5) 的 0.9 倍系数让步长自然增长。4. 容限、初始步长和缩放因子边界设错了照样跑不动4.1 绝对容限和相对容限怎么配合自适应变步长误差阈值 tol 在很多教材里只有一个数字但落到 MATLAB 里要区分绝对容限和相对容限。绝对容限 atol 决定解接近 0 时允许的绝对偏差相对容限 rtol 决定解量级较大时允许的相对偏差。例如在电路仿真里电压和电流两个状态量量级差上百倍如果只用单一 abs 误差电流分量会被误判成误差很大程序会不断把步长缩小最后卡死在初始化阶段。常见的处理是沿用 ode45 的参数习惯rtol 1e-6; atol 1e-8; err sqrt(sum(((y_half - y_full) ./ (atol rtol*abs(y_half))).^2));注意 atol 和 rtol 都是一个标量作用于所有分量这其实是一种简化。更细的用法是让 atol 和 rtol 与 y 同维度对不同分量分别设置容限。比如初值问题中某个状态只变化 1e-5 量级那就单独给它一个 atol1e-7而另一个状态变化到 1e3 量级仍用 rtol 控制。分量级容限能显著减少刚性方程组里的无效拒绝。4.2 初始步长 h0 的经验做法自适应步长理论上是自动找步长但初始步长给得离谱会让循环在开始阶段空转几轮。最稳妥的经验是先跑一步欧拉试探取一个很小的临时步长 h_try算出 dydt 的初值用 dydt 变化率反推初始 h。对罗森布罗克这类刚性求解器更常见的策略是 h0 sqrt(eps) * max(norm(y0), 1e-6)。对本项目这种基于 RK4 的代码我一般直接给h0 1e-3 * max(norm(y0), 1);这个值对大多数非病态问题能在前几步内快速放大到合理区间。如果第一步就反复拒绝则说明 h0 比可接受范围大太多此时不要手动把 h0 调小到离谱的程度而是先检查 f.m 里的时间尺度。比如方程系数是 1e6而 h0 仍从 1e-3 出发第一步 f(t,y) 就会让 k2 和 k4 相差好几个数量级误差信号失真后面调整等于在噪声上做控制。更合理的做法是按 f 中最大时间常数的量级给初始步长。4.3 缩放系数、最小步长和最大步长的钳制步长更新公式 min(max(0.9 * (tol/err)^(1/5), low), high) 里的 low 和 high 不只是保护误差也保护程序本身。low 设置太小例如 1e-10一旦误差信号有异常峰值h 直接缩到接近机器精度后面就算误差恢复正常也要几百步才爬回来。更常见的做法是用 H 序列里最小步长做一个全局钳制当 h 1e-8 * tspan 时直接触发“求解失败”退出而不是继续死循环。if h 1e-8 * abs(tspan(2) - tspan(1)) warning(步长低于时间尺度的1e-8可能已经遇到奇点); break; end这里 1e-8 不是需要反复调节的神秘系数而是基于双精度浮点余量的保守选择。如果方程本身是刚性最小步长会在合法区间内自动降下来如果 h 降到这个阈值以下说明右端函数在这点附近有奇点或非连续再自适应也没有意义。high 上限方面不赞成把 high 设成 10 以上。自适应步长通过误差估计放大步长但如果放大速度过快区间内快速变化会被跨过去。比如 h 从 1e-4 放大到 1e-3放大 10 倍后误差也许刚好在阈值内但下一步误差突然超限程序又要缩回来形成一种低频振荡。上限 5.0 配合 0.9 的安全系数让 h 的增长路径更平滑。下面这张表可以作为参数调整入口调整目标改哪个参数副作用注意让结果更精确减小 atol / rtol总步数线性上涨刚性区间可能涨几倍让平缓段更快增大 high可能在快慢过渡区错过一个尖峰减少拒绝次数增大安全系数到 0.95单步误差变大但总调用次数可能下降解决启动失败改 h0 数量级只影响前几步不会改变后续步长轨迹最后要提醒在 zhukong.m 里错误地将最大步长 high 设成 100 并不会让程序跑得更快它只会让 t 推进的跳变变大然后每次误差检查都会失败最终总计算量反而增加。自适应方法的核心价值在于“该小的地方小该大的地方大”maxstep 上限只是一个防呆边界。5. 验算自适应行为用 H 输出序列反推刚区间位置5.1 把调参结果可视化zhukong.m 返回的 H 数组就是每步实际使用的步长。运行完求解后画两条线一条是解 y(t)一条是 h(t)刚区间立刻显形[T, Y, H] zhukong(f, [0 5], 1, 2e-4, 1e-6); subplot(2,1,1); plot(T, Y); title(y(t)); grid on; subplot(2,1,2); semilogy(T(2:end), H(2:end)); title(step size h(t)); grid on;这段代码用 semilogy 而不是 plot因为步长在最初几十步可能横跨 1e-5 到 1e-1对数纵轴才能看出指数级变化。如果发现 h(t) 在某个区间出现锯齿形骤降骤升说明误差在阈值附近来回震荡优先检查 err 那行的分母有没有用 atol rtol*abs(y_half)。5.2 与 ode45 做同点差值验证自适应实现写完后最直接的正确性验证是插值后和 ode45 结果对照。ode45 内部使用带误差控制的 Dormand-Prince 方法这套“全步长减半步长”是一对更朴素的策略两者不应该完全相同但应具有相近的数量级[T2, Y2] ode45(f, [0 5], 1, odeset(RelTol, 1e-6, AbsTol, 1e-8)); y_interp interp1(T, Y, T2, pchip); max(abs(y_interp - Y2))max 差值如果在 1e-4 量级以下说明主循环没有原则性错误。如果差值达到 1e-2第一步不是调小容限而是看 H 序列有没有大跳跃尤其关注 t0 附近前 10 步。多数实现问题都出在拒绝分支里忘记把 h 更新回下次迭代使用的变量导致实际步长和 H 记录不一致。5.3 把“拒绝次数”变成调参指标再给一个能直接量化的技巧在 zhukong.m 的 else 分支里加一个计数器统计整个求解过程中步长被拒绝的次数。这个数字比最终步数更能说明问题。容限设置合理时拒绝次数应远小于接受次数如果拒绝次数超过接受次数的 20%说明初始步长偏大或安全系数过小应当把 0.9 提到 0.95而不是把 tol 放大。相反如果 H 序列几乎不变说明步长一直被上限钳住此时应当调大 high或确认方程是否真的存在刚性区间。这一行放在 else 分支末尾即可reject_count reject_count 1;在函数返回前把 reject_count 作为第三个输出或写入全局变量。用静态分析看一个积分算法是好是坏单看总步数并不准确用“接受步数 / 总评估次数”来衡量更接近真实成本。half.m 每次调用会执行四次 f 求值一次全步长加两次半步长总计 12 次求值这个数才是和 ode45 对比性能时应关注的核心数值。本文还有配套的精品资源点击获取
