拆解ode45：Runge-Kutta四阶定步长实现与误差分析

简介：围绕Matlab中龙格-库塔（Runge-Kutta）方法展开的这份文档，面向需要求解常微分方程数值解的学生、科研人员与工程实践者，重点梳理ode45所采用的四阶五阶Runge-Kutta单步算法原理，并给出自行编写定步长Runge-Kutta函数的实现思路。压缩包共1个doc文件，体积约70KB，内容以算法说明、ode45调用方式、一阶与高阶方程示例及自编函数对照为主，适合作为数值分析课程学习与Matlab编程的补充材料。已有4541人学习下载，说明其在常微分方程数值解入门中具有一定参考价值。读者可从中了解k1至k4的中间斜率计算、加权平均更新公式、ode45参数与options设置、刚性问题换用ode23的判断，以及高阶方程转一阶方程组的处理方式，并对照自编runge_kutta1函数理解底层实现与内置函数之间的差异。

1. 一行 ode45 背后藏着的公式：这套 Runge-Kutta 资源到底补的是什么

在 Matlab 里解常微分方程，多数人第一反应就是敲一行[t,y] = ode45(...)，出图就算收工。可一旦曲线跟手算对不上、步长改小结果反而变了、或者命令行十几秒没有输出，这个黑匣子就开始咬人。Runge-Kutta 方法就是压在 ode45 底下的那套公式，ode45 用的是四阶五阶嵌入式系数对，靠两组结果之差估误差，再回头决定这一步走多大。这份资源补的正是被封装掉的那一层：k1 到 k4 怎么算、1/6 的权重从哪来、定步长版本自己怎么写、高阶方程怎么降阶、自编结果怎么跟 ode45 对齐。适合已经会调 ode45 但说不清它在干什么的人，也适合要把同一套积分搬到 C 端或嵌入式里、需要一个可逐行对照实现的工程师。

2. 拆开 k1 到 k4：四阶 Runge-Kutta 的加权斜率与定步长实现

2.1 从泰勒展开到平均斜率：四个采样点为什么是这四个

一阶方程写作y' = f(x,y)，最朴素的离散化是前向欧拉：把导数当成常数，y(n+1) = y(n) + h*f(x(n), y(n))。这一步只用了区间左端一个点的斜率，局部截断误差是 O(h²)，跑长区间时会肉眼可见地往下掉。改进欧拉（也就是二阶 Runge-Kutta）在区间中点再采一次斜率，用两次斜率的平均去推下一步，局部误差降到 O(h³)。四阶的思路一样，只是采样点从两个扩到四个：左端点、两个中点（用不同方式预测）、右端点，然后按 1:2:2:1 加权。

从微分中值定理看会更清楚：y(n+1) = y(n) + h*K，其中K = f(x(n)+θh, y(x(n)+θh))是区间内某个未知点的斜率，0<θ<1。整个 Runge-Kutta 家族干的事，就是用若干个已知点的 f 值去凑这个 K。二阶凑一次中点，四阶凑四次，通过选择系数让展开式与泰勒级数在 h⁴ 项之前完全吻合，于是局部截断误差落到 O(h⁵)，全局误差 O(h⁴)。资源里给出的那组经典系数就是：

K1 = f(Xn, Yn) K2 = f(Xn + h/2, Yn + (h/2)*K1) K3 = f(Xn + h/2, Yn + (h/2)*K2) K4 = f(Xn + h, Yn + h*K3) Yn+1 = Yn + h*(K1 + 2*K2 + 2*K3 + K4)/6

注意 K2 和 K3 都在半步处取值，区别在于预测时用的斜率不同：K2 用左端点的 K1 外推，K3 用 K2 再修一次。这就是所谓「预估—校正」在单步内的体现，也是四阶比二阶贵一倍的代价所在——每前进一步要算四次右端函数。

这里有个容易被忽略的边界：四阶 RK 的精度建立在 f 足够光滑的前提上。如果右端函数在某些点有跳跃、绝对值、分段定义，泰勒展开的前提就破了，理论上的四阶在实际里可能连二阶都不如。我遇到过带符号函数的摩擦模型，ode45 和自己写的定步长 RK 都在切换点附近抖，最后是把切换点当成事件、分段积分才稳住。这一点在选方法时比阶数更重要。

2.2 手写 runge_kutta1：签名、步数计算与列向量预分配

资源里给的定步长实现只有二十来行，但每一行都有讲究。我按原样复现一遍，顺手补上预分配和维度强制，避免跑大区间时被内存重排拖慢：

function [x, y] = runge_kutta1(ufunc, y0, h, a, b) % ufunc : 微分方程组右端函数句柄，签名 dydt = ufunc(t, y) % y0 : 初值列向量，n×1 % h : 定步长 % a, b : 积分区间起点与终点 n = floor((b-a)/h); % 完整步数，末尾不足一步的部分直接舍掉 x = zeros(1, n+1); % 时间序列预分配 y = zeros(numel(y0), n+1); % 解矩阵预分配，行数 = 状态维数 x(1) = a; y(:,1) = y0(:); % 冒号强制成列向量，防止行向量传进来维度错乱 for ii = 1:n x(ii+1) = x(ii) + h; k1 = ufunc(x(ii), y(:,ii)); k2 = ufunc(x(ii)+h/2, y(:,ii) + h*k1/2); k3 = ufunc(x(ii)+h/2, y(:,ii) + h*k2/2); k4 = ufunc(x(ii)+h, y(:,ii) + h*k3); y(:,ii+1) = y(:,ii) + h*(k1 + 2*k2 + 2*k3 + k4)/6; end end

逻辑上有三处值得说。第一，n = floor((b-a)/h)是向下取整，意味着积分终点不一定正好落在 b 上，比如(4-0)/0.3得 13.33，实际只走到 3.9。要做等间距对拍，就把(b-a)/h设成整数，或者改成linspace(a, b, n+1)并让 h 由(b-a)/n反算。第二，y0(:)那一下看着多余，但当调用方从别的函数里传进来一个行向量时，y(:,ii)会静默产生一个 1×n 的行，下一轮h*k1/2的广播规则就会把它变成矩阵，错误要到几十步之后才炸，非常难查。第三，四个 k 都是一次函数调用，ufunc里最好不要有persistent变量或打印语句，否则每步四次调用会把开销放大到不可接受。

参数上，h是唯一能控制精度的旋钮，但它不是越小越好。四阶 RK 的全局误差理论上是 O(h⁴)，h 减半误差应该掉到 1/16 左右；可当 h 小到 1e-6 量级，双精度的舍入误差开始累积，步数上百万之后误差反而回升。实际做工程，我一般先用h = (b-a)/1000起步，看曲线形状对不对，再减半做一次收敛验证，确认误差确实按 16 倍在降，才敢把结果拿去用。

2.3 三个参数怎么定：h、区间端点与初值维数

h的选取跟系统的特征时间尺度挂钩，不是拍脑袋。粗略的判据是h要小于最快变化周期的一个零头，经验上取特征时间的五十分之一到百分之一，同时保证一个周期内至少采到二十个点，否则画出来的曲线会假性光滑。资源里的例子用h=0.25跑[0,15]，60 步出结果，曲线形状没问题，但如果要跟 ode45 做逐点误差对比，0.25 的粗步长已经能看出几百个微单位的偏差了，得降到 0.01 量级才谈得上「对齐」。

区间端点a、b决定步数，也决定后面画图的横轴。有一点常被忽略：定步长 RK 输出的时间点是自己累加出来的，而 ode45 输出的 T 是自适应步长下的非均匀点，两边的采样位置根本对不上。直接拿plot(T,F)和plot(T1,F1)并排看还行，要做norm(F - F1)这种数值对比，必须先用interp1把其中一条插到另一条的网格上，或者干脆指定odeset('MaxStep', h)让 ode45 的步长上限不超过 h，输出点会更密。

初值维数是最容易翻车的地方。y0的维数必须严格等于右端函数返回的列向量长度，二者不一致时 Matlab 不会报「维数不匹配」，而是报Unable to perform assignment because the size of the left side is 2-by-1 and the size of the right side is 3-by-1，或者更糟——在广播规则下静默算出一个错误的结果。资源里特别提了一句「测试时改变 test_fun 的函数维数，别忘记改变初始值的维数」，就是踩过这个坑之后写的。

3. 高阶方程降阶：把 y'' 写成状态向量再交给 RK

3.1 降阶的机械操作：y1=y，y2=y'，右端写成列向量

所有一阶以上的常微分方程，进 RK 之前都必须先降成一阶方程组，这一步没有例外。做法固定：设y1 = y，y2 = y'，y3 = y''，一直到yn = y(n-1)，那么原来那个 n 阶方程F(y, y', ..., y(n), t) = 0就被拆成 n 个一阶方程，第一个是y1' = y2，第二个是y2' = y3，最后一个由原方程解出最高阶导数得到。写成矩阵就是dydt = [y2; y3; ...; f(t, y1, y2, ..., yn)]。

这个操作本身机械，出错集中在两处。一处是最后一个分量的符号：原方程里的y''系数是 2 的时候，解出来的y''要除以 2，漏掉这个除法是最常见的错误，而且结果看起来「差不多」，只有做误差对比才暴露。另一处是初值的对应关系：y0 = [y(0); y'(0)]，顺序不能反，反了之后得到的解可能同样光滑、同样收敛，但完全是另一条曲线。

3.2 一个二阶例子：从 y'' = -t·y + e^t·y' + 3sin(2t) 到 odefun

拿资源里那个二阶方程落地。原式y'' = -t*y + e^t*y' + 3*sin(2t)，令x(1) = y、x(2) = y'，右端函数写成：

function dx = test_fun(t, x) % 二阶方程 y'' = -t*y + exp(t)*y' + 3*sin(2t) % 状态定义：x(1) = y，x(2) = y' dx = zeros(2,1); % 先定列向量，别用 [] 在循环里拼 dx(1) = x(2); % y' = x(2) dx(2) = -t*x(1) + exp(t)*x(2) + 3*sin(2*t); % 由原方程解出 y'' end

两个细节。dx = zeros(2,1)必须先写，因为后面dx(1)的赋值在空矩阵上会创建一个 1×1 的标量，接着dx(2)赋值时 Matlab 补零扩展成 1×2 的行向量，返回维度就是错的，RK 那边h*k1/2一乘，广播出 2×2 矩阵，报错信息完全指不到源头。另外exp(t)在 t 较大时会迅速吃掉精度，积分区间拉到 10 以上时建议先做变量替换或缩放时间，否则 ode45 会把步长压得很小、速度骤降。

调用侧，ode45 和自编 RK 用的是同一个函数句柄，这一点很重要——只要右端函数是同一个，两条曲线才有可比性：

tspan = [0 4]; % 资源原例写的是 [3.9 4.0]，区间只有 0.1，曲线近似直线 y0 = [2; 8]; % 对应 y(0)=2, y'(0)=8，顺序不能反 [T, X] = ode45(@test_fun, tspan, y0); [T1, X1] = runge_kutta1(@test_fun, y0, 0.01, tspan(1), tspan(2)); subplot(1,2,1); plot(T, X(:,1), '-o'); title('ode45'); legend('y','y'''); subplot(1,2,2); plot(T1, X1(1,:), '-*'); title('自编 RK4');

X(:,1)取的是状态矩阵的第一列，也就是 y 本身，第二列是 y'。资源原例里用的是plot(t, x(:,1), '-o', t, x(:,2), '-*')，一次把两个状态画出来，好处是能同时看到位移和速度是否自洽——如果 y 光滑而 y' 出现锯齿，八成是右端函数写错了。

3.3 自编 RK 与 ode45 并排跑：步长、采样点与误差量级

两个函数都跑完之后，别只看图。图会把 1e-3 量级的偏差抹平，看着「两条曲线重合」就以为对了，实际可能差着一个数量级。我的习惯是固定终点，直接算末值差：

yEnd_ode = X(end, :); % ode45 末值（1×2） yEnd_rk = X1(:, end).'; % 自编 RK 末值，转成 1×2 方便相减 relErr = norm(yEnd_rk - yEnd_ode) / norm(yEnd_ode); fprintf('末值相对误差 = %.3e\n', relErr);

h = 0.01时，四阶 RK 在这个例子上的末值相对误差通常落在 1e-8 到 1e-6 之间；把 h 改成 0.001，误差应该再降三到四个数量级。如果降不下来，先怀疑三件事：初值顺序、右端函数的符号、y0的维度。这三样对不上时，误差不会发散，而是稳定地停在一个百分位，非常像「算法精度不够」，很容易误判成 RK 本身的问题。

另外提醒一句，ode45 默认RelTol = 1e-3，它自己的解也不是「真值」。拿它当基准做对比时，要么把odeset('RelTol',1e-10,'AbsTol',1e-12)调紧，要么用h极小的自编 RK 当参考解，否则你测的是「两种误差的差」，不是算法误差。

4. ode45 的 options 与事件：RelTol、AbsTol 和 sol 结构体怎么用

4.1 odeset 里真正该动的三个参数

[T,Y] = ode45(odefun, tspan, y0)是最常见的形式，但它丢掉了 ode45 最有价值的部分——求解控制和事件检测。完整签名是[T,Y] = ode45(odefun, tspan, y0, options)，options由odeset构造。参数表里能设的项有几十个，实际工程里我会动的只有三个：

写法是options = odeset('RelTol',1e-8,'AbsTol',1e-10,'MaxStep',0.01);，然后传给 ode45。RelTol调太紧的直接后果是步长被压到极小、步数暴涨，一个本来两秒跑完的模型能卡上几分钟；而AbsTol设得比解的幅值还大时，误差判据形同虚设，求解器会大摇大摆地跨过细节。我一般先用默认值跑通，确认曲线形态没问题，再按需收紧一档，对比两次结果差异，差异可接受就停手。

4.2 事件函数：过零检测怎么写才不白写

带事件的求解是这个资源里最实用、也最少人用的部分。典型场景是撞地、开关切换、达到阈值停机——你只关心「第一次发生」的时刻，不想让求解器一路算到底再回头找。事件函数签名固定为三返回值：

function [value, isterminal, direction] = evt_ground(t, y) % y(1) 代表高度，检测它从正穿到负的瞬间 value = y(1); % 过零即触发 isterminal = 1; % 1 = 触发后立即停止积分；0 = 只记录不停止 direction = -1; % -1 = 只捕捉由正到负的穿越；+1 反之；0 = 两个方向都要 end

挂上去的写法是options = odeset('Events', @evt_ground);，然后[T,Y,TE,YE,IE] = ode45(@odefun, tspan, y0, options);。TE 是事件时刻，YE 是对应状态，IE 是事件编号——多个事件函数打包成元胞数组传入时，IE 告诉你触发的是第几个。资源正文里把 TE/YE/IE 写成了「事件发生时间 / 事件解决时间 / 事件消失时间」，后两个说法容易被误读，按 Matlab 官方语义，YE 是事件发生时刻的状态值，IE 是事件索引。

direction是最容易写错的一项。设成 0 表示两个方向都捕捉，弹跳问题里小球落地又弹起，会在同一个物理位置反复触发，如果isterminal = 1，积分在第一次落地就停了，后面的弹跳全丢；正确做法是isterminal = 0加direction = -1，让它只记落地时刻，然后手动改初值重启下一段。

4.3 sol 结构体：不等间距输出怎么取任意时刻的解

sol = ode45(odefun, [t0 tf], y0)返回的是结构体而不是矩阵，里面存的是求解器内部的步点和插值系数。好处是不用逼求解器在你关心的每个时刻都落点，取解时用deval：

sol = ode45(@test_fun, [0 4], [2; 8], odeset('RelTol',1e-8)); tq = linspace(0, 4, 200); % 想在哪看就在哪看，与内部步长无关 Yq = deval(sol, tq); % 2×200，每列是一个时刻的状态 plot(tq, Yq(1,:), tq, Yq(2,:));

这比[T,Y]强的点在于：输出网格完全由你定，不会因为求解器在某个陡峭段把步长压到 1e-6 而拿到一堆挤在一起的点；同时deval用的是求解器自带的稠密输出插值，精度是四阶的，比事后interp1的线性插值靠谱得多。要注意deval只在[t0, tf]内有定义，超出范围会报错，需要外推就得自己接。

5. 自编 RK 和 ode45 对不上：五类排查现场

5.1 曲线形状对但幅值差一截：初值维数和行列向量

现象：两条曲线趋势一致，峰值位置也对，但整体差百分之几，改步长也不收敛。

原因：y0传的是行向量（比如[1 1 3]），进入runge_kutta1后y(:,1) = y0(:)虽然强制成列，但调用方在别处仍然按行向量理解，或者在plot(T1, F1)里把状态矩阵转置了。另一种是y0的维数比右端函数返回的少一项，Matlab 在某些版本下用广播补零，不报错。

解决：所有初值一律写成列向量[1; 1; 3]，函数内部第一行打印size(y0)自查；状态矩阵统一约定「行是维数、列是时间」，画图前用size(X1)确认是n×N而不是N×n。

5.2 自编结果一开始准、后面越跑越偏：定步长与累积误差

现象：前几十步和 ode45 几乎重合，跑到后半段开始肉眼可见地分开，且 h 越小分开得越晚。

原因：定步长 RK 的每一步误差都独立叠加，全局误差按 O(h⁴) 累积；而 ode45 是自适应步长，误差超阈值就自动缩步，所以长区间上天然占优。这不是自编代码写错了，是方法本身的边界。

解决：做算法验证时把区间截短（比如[0 2]而不是[0 15]），或者把 h 降到(b-a)/5000再看；做长期仿真就老实换 ode45 或 ode15s，别指望定步长 RK 扛住长区间。

5.3 跑很久没输出或报积分容差失败：刚性问题

现象：ode45 卡住不出结果，或者报Unable to meet integration tolerances without reducing the step size below the smallest value allowed。

原因：方程是刚性的——解里同时存在快变和慢变分量，显式 RK 的稳定域要求步长小于快变尺度，可总区间又由慢变尺度决定，于是步数爆炸。资源正文里说「换用 ode23 来解」，这条建议我认为偏保守：ode23 是低阶显式方法，稳定域比 ode45 稍大但仍然是显式，对真正的刚性问题一样会卡，只是卡得晚一点。

解决：换ode15s（变阶 BDF/NDF，刚性问题的主力）或ode23s（二阶 Rosenbrock，适合中等刚性且精度要求不高），配上解析或数值 Jacobian 会快很多。判断刚性有个土办法：把右端函数的雅可比矩阵在初值附近求特征值，最大最小模之比超过两三个数量级，基本就是刚性。

5.4 报下标越界：odefun 里 y 的索引写错

现象：Index exceeds the number of array elements，或者Unable to perform assignment because the left and right sides have a different number of elements。

原因：降阶时状态编号和方程编号错位，比如三阶方程该有y(1)、y(2)、y(3)，代码里写成了y(0)（Matlab 下标从 1 开始），或者dy = zeros(3,1)之后只赋了前两个分量，第三个留零——后者不报错，但解是错的。

解决：右端函数第一行写assert(numel(y) == 3, '状态维数应为 3')，把错误挡在入口；赋值顺序严格按dy(1)=y(2); dy(2)=y(3); dy(3)=...写，不要跳号。

5.5 事件函数写好了却一次都不触发

现象：明明曲线穿过了阈值，TE返回空，或者IE只有第一个事件。

原因：direction符号反了，或者value写成了恒正量（比如忘了减去阈值），或者多个事件函数用元胞数组传入时isterminal全设成 1，第一个事件一触发就整体停机。

解决：把value单独画出来看一眼过零点；direction先设 0 排除方向问题，确认能触发后再收窄；多事件场景下只给需要停机的那个设isterminal = 1，其余设 0。

6. 收敛阶验证与刚性切换：两个立刻能做的动作

6.1 用减半步长法把「四阶」验出来

四阶 RK 的全局误差是 O(h⁴)，这意味着 h 减半，误差应该掉到原来的 1/16，取以 2 为底的对数就是 4。这个检验只要五行代码，比看曲线可靠得多：

hlist = [0.2 0.1 0.05 0.025 0.0125]; err = zeros(size(hlist)); for k = 1:numel(hlist) [~, Y] = runge_kutta1(@test_fun, [2; 8], hlist(k), 0, 1); err(k) = abs(Y(1,end) - yRef); % yRef 用极小步长或紧容差 ode45 算出的末值 end p = log2(err(1:end-1) ./ err(2:end)); % 应稳定在 4 附近 disp(p)

yRef一定要用比最小 h 还精细的解，否则测的是参考解自身的误差。我一般用h = 1e-4的自编 RK 或odeset('RelTol',1e-12)的 ode45 当基准。如果 p 算出来在 2 附近，说明右端函数不光滑或者实现里某一步系数写错了；如果 p 随 h 减小反而下降，那是舍入误差开始主导，说明 h 已经太小，该停手了。这个动作我每次改完 RK 实现都跑一遍，比对着公式一行行核系数快得多。

6.2 刚性识别与求解器切换的对照表

把资源里提到的几个 ode 函数放一张表里对齐，选型时直接查：

切换的动作很简单：把ode45换成ode15s，其他签名一字不改。真正要花时间的是判断到底该不该换——我的判断顺序是先看 ode45 是否报容差失败，再看解里有没有相差三个数量级以上的时间尺度，最后才动特征值。这三步都没问题就别乱换，隐式方法每步要解线性方程组，小模型上反而更慢。

有一回我图省事，把一个明摆着刚性的化学反应模型丢给 ode45，跑了一晚上没出结果，第二天换成 ode15s 三十秒结束，白等的那一宿就是没做特征值估计的学费。从那以后我每次拿到一组新方程，都先用初值附近的雅可比估一遍特征值量级，再决定用显式还是隐式，不再凭感觉试。希望帮到你。

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