简介：面向MATLAB数值计算初学者与需要求解常微分方程的工程学习者，这份文档围绕龙格-库塔（Runge-Kutta）方法展开，重点说明ode45所采用的四阶五阶单步算法、截断误差特点及非刚性问题的求解思路。压缩包共1个doc文件，大小约70KB，内容紧凑，便于作为算法入门、课程复习与代码查阅的便携资料。文档梳理了ode45的调用形式、参数含义与odeset选项设置，并结合一阶、高阶常微分方程给出变量替换、求解与绘图示例；同时提供自编定步长龙格-库塔函数的完整实现，便于对照ode45理解k1至k4的计算流程、加权更新公式及与自带函数的差异。内容还涉及事件检测输出、相关ode函数的选择，以及对刚性与非刚性问题的判别提示。已有4541人学习，适合希望从数学原理到MATLAB落地逐步对照的读者。
1. 手写 RK4 和 ode45 摆在一起：先看清两条曲线差在哪
三方程耦合的常微分方程组，初值[1 1 3]，时间区间[0 15]。同一份右端函数，一半交给 Matlab 自带的ode45，一半喂给自己写的四阶 Runge-Kutta 定步长循环，两条曲线画在同一张图上，形状大体贴合，但采样点数量差了十几倍，末端还悄悄分叉。很多人对数值解 ODE 的第一印象就卡在这个画面上：ode45凭什么敢用变步长，手写的 Runge-Kutta 为什么只能老老实实一步步挪。
这份资料把两件事摊开讲了：四阶 Runge-Kutta 的四个斜率怎么来的、ode系列函数的调用签名每一项填什么。它适合正在用 Matlab 解微分方程、却被步长、容差、刚性这些词反复绊住的人，也适合想把ode45当黑匣子用、结果一出问题就只能干瞪眼的人。下面按“公式怎么落地 → 接口怎么调 → 高阶怎么降 → 哪里会翻车”走一遍。
2. 四阶 Runge-Kutta 的四次斜率：把公式拆成能跑的循环
2.1 从欧拉的一步到 RK4 的四步
一阶常微分方程写成y' = f(x,y)，最朴素的做法是把导数当成差商：(Y(n+1) - Yn)/h ≈ f(Xn,Yn)，移项就得到前向欧拉公式Y(n+1) = Yn + h*f(Xn,Yn)。这个式子的几何含义很直白：拿当前点的斜率，笔直往前走一步。问题也直白，它只用了区间左端一个点的信息，相当于用矩形面积代替曲边梯形，局部误差是O(h²)，步长稍微放大，曲线立刻塌。
微分中值定理给了另一条路：存在某个0 < t < 1，使得Y(n+1) = Yn + h*f(Xn + t*h, Y(Xn + t*h))。这里的f值叫区间上的平均斜率，欧拉法粗暴地取了t = 0。Runge-Kutta 的思路是：既然一个点估不准平均斜率，那就在区间内多预计算几个点，把它们的斜率加权平均。取两个点就是二阶，取四个点就是四阶。
四阶的四个斜率是这么排布的：
| 斜率 | 计算位置 | 依赖 |
|---|---|---|
| K1 | 区间左端点(Xn, Yn) | 当前解 |
| K2 | 半步长点，用 K1 推进 | K1 |
| K3 | 半步长点，用 K2 推进 | K2 |
| K4 | 整步长点，用 K3 推进 | K3 |
合成公式是Y(n+1) = Yn + h*(K1 + 2*K2 + 2*K3 + K4)/6。权重1/6、2/6、2/6、1/6不是拍脑袋凑的，是把Y(n+1)在Xn处做泰勒展开，逐项和f的各阶偏导匹配出来的结果，匹配到h⁴项，局部截断误差就落到O(h⁵)，全局误差O(h⁴)。中间那两步推导确实繁琐，工程上记住两条就够用：每步四次函数调用，换四阶精度；局部误差是步长的五次方，步长减半，单步误差降到 1/32。
顺带说一句原料里的表述：“ode45采用四阶五阶 Runge-Kutta 单步算法，截断误差为(Δx)³”。这个三次方的说法是简化，ode45实际用的是 Dormand-Prince 4(5) 对，用五阶公式估计误差、用四阶公式推进，工程上按四阶全局精度理解不会出错，但别拿三次方去反推步长。
2.2 定步长 runge_kutta1：参数表、预分配与列向量约定
把上面的公式写成函数不难，难的是写成能同时吃标量和向量的版本。资料里给出的runge_kutta1就是这种写法，我按自己的习惯补了预分配和维度保护：
function [x, y] = runge_kutta1(ufunc, y0, h, a, b) % ufunc 微分方程右端函数句柄，签名 dy = ufunc(x, y)，dy 必须是列向量 % y0 初始值，标量或列向量，长度等于方程个数 % h 定步长，正数 % a, b 时间起点与终点，要求 b > a 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 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; x(ii+1) = x(ii) + h; end end逐项说清楚。n = floor((b-a)/h)意味着(b-a)不是h的整数倍时，最后那一小段不会被积到，想要严格落在b上，就得自己把h调成能整除的长度，或者把最后一步单独处理。zeros(numel(y0), n+1)这一行不能省，Matlab 的数组在循环里逐列追加会反复重新分配内存，几百步看不出来，几十万步就是肉眼可见的慢。
y0(:)这个写法是血泪经验。调用方传[1 1 3]进来是 1×3 的行向量，如果不转置，y(:,1) = y0在部分版本里会报维度错误，在另一些版本里会悄悄广播成 3×3，后面所有k值全错，图还能画出来，只是完全不对。函数内部统一列向量，调用方就不用每次记得加转置。
拿资料里的一阶方程试一下：
odefun = @(t, y) (y + 3*t) / t^2; % 方程 t^2*y' = y + 3t [t, y] = runge_kutta1(odefun, -2, 0.01, 1, 4); plot(t, y, '-'); xlabel('t'); ylabel('y'); title('t^2y'' = y + 3t, y(1) = -2');这里y0 = -2是标量，numel得 1，y是 1×(n+1) 的行向量，y(:,ii)取出的是标量，k1到k4也都是标量，整条链路退化成纯量运算，跑得通。要验证精度，拿解析解(3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)在同一批t上求值做差即可，Ei是指数积分函数，Matlab 里写作ei。步长 0.01、区间长度 3，共 300 步，末端误差量级通常在1e-8上下，换个h = 0.02再跑一遍，误差应该涨到十几倍，这个比值就是四阶收敛最直观的证据。
2.3 步长 h 不是越小越好：截断误差与舍入误差的拉锯
新手最容易形成的直觉是“步长越小越准”，这句话只对了一半。定步长 RK4 的总误差由两部分叠加：截断误差随h减小按h⁴下降，舍入误差随步数增加按1/h上升，两条曲线一降一升，中间有个最优点。步长小到1e-8以下，双精度浮点的有效位被逐步蚕食，结果反而变差，曲线开始出现细碎的毛刺，这就是典型的“算得越认真越不准”。
实践经验：先按(b-a)/100到(b-a)/1000取一个h，跑一遍看曲线形状；再做步长减半对照，如果末端解的相对变化已经小于你能接受的精度，就停手。别一上来就1e-6，既慢又未必更准。
还有一个更隐蔽的约束是稳定性，跟精度无关。显式 RK4 的绝对稳定域在实轴负方向大致到-2.78附近，意思是方程里出现快速衰减模态、且λ*h的模超过这个量级时，数值解会震荡发散，哪怕真解光滑得很。这时候把h再调小反而“治好”了，其实是把λ*h重新压回稳定域里，代价是步数暴涨。这类方程就是后面要说的刚性方程，正确做法是换隐式求解器，不是硬压步长。判断方式很简单：把h减半，如果误差按 16 倍降，那是精度主导；如果误差纹丝不动甚至变大，而曲线还在抖，那就是撞上稳定性边界了。
3. ode45 的三个入参和一堆 options：调用签名逐项拆
3.1 odefun、tspan、y0 的合法写法和错位重灾区
最常用的调用形式是[T, Y] = ode45(odefun, tspan, y0)，三个入参各有各的坑。
odefun是右端函数的句柄，三种写法都合法：函数文件名要加@，比如@test_fun；匿名函数直接写@(t,y) ...；老版本的内联函数inline(...)现在基本不用了，能跑但性能差。注意@test_fun和test_fun是两回事，后者会把函数当变量求值，报错信息还特别含糊。函数签名固定是dy = odefun(t, y)，参数顺序不能反，t在前y在后，跟数学写法f(x,y)的直觉正好相反，这是最高频的错位点。
tspan有两种形态，行为完全不同。写成两元素[t0 tf]时，求解器按自己的自适应步长内部推进，返回它实际踩过的所有节点，点数由求解器决定，可能几十也可能几千。写成多元素向量[t0, t1, ..., tf]时，求解器仍然按内部步长算，但只在指定的这些时刻返回插值结果，点数是固定的。想要等间距输出又不想干预内部步长，用第二种；想让求解器自己决定采样密度，用第一种。
odefun = @(t, y) (y + 3*t) / t^2; tspan = [1 4]; y0 = -2; [T, Y] = ode45(odefun, tspan, y0); % 自适应节点，T 是列向量 plot(T, Y, '-o');y0必须是向量，长度等于方程个数，标量方程就传标量。传行向量[1 1 3]在多数版本里能跑，但内部会做一次形状转换，一旦你的odefun返回的也是行向量，两者对不上就会报维度错误。统一用列向量[1; 1; 3]是最省心的写法。
返回的T是列向量，长度和Y的行数一致；Y是矩阵，每一列对应一个分量，取第一个分量写Y(:,1)。这一点跟plot(T,Y)直接画多条线是自洽的，但和很多人的“每行一个时刻”的直觉相反，取错维度画出来的图会是一堆乱七八糟的折线。
3.2 odeset 里真正影响结果的四个参数
options用odeset构造，参数名不区分大小写，但拼错不会报错，只会被静默忽略，这是个经典的坑。真正常用的就那么几个：
| 参数 | 默认值 | 作用 | 什么时候改 |
|---|---|---|---|
RelTol | 1e-3 | 相对误差容限 | 要更高精度时降到1e-6到1e-10 |
AbsTol | 1e-6 | 绝对误差容限，可标量或向量 | 分量量级差异大时按分量给 |
MaxStep | 区间长度的 1/10 | 步长上限 | 怕错过快速变化段时压小 |
Events | 空 | 事件函数句柄 | 需要抓过零、终止积分时 |
误差判据是“或”的关系：第i个分量的估计误差满足err_i <= max(RelTol*|y_i|, AbsTol_i)就算通过，任一条件满足即可。这意味着当某个分量数值接近零时，RelTol那一项自动失效，全靠AbsTol兜底。如果你的解里有个分量在1e-8量级附近晃，默认AbsTol = 1e-6等于完全不管它，曲线会呈现明显的阶梯状，很多人以为是方程有问题，其实是容差放得太松。
MaxStep的作用容易被低估。求解器在解变化平缓的区域会主动放大步长，理论上误差控得住，但如果方程里有窄脉冲、或者你后续要拿结果做频谱分析，采样点太稀会直接丢掉高频成分。把它压到0.01这种量级，代价是步数变多，收益是曲线细节保得住。
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.05); [T, Y] = ode45(odefun, [1 4], -2, opts);改容差要看收益比。RelTol从1e-3降到1e-6，步数大概翻几倍；再降到1e-10，步数继续涨，但结果的可见精度提升已经很有限，除非你在做长期积分或者对初值敏感的问题，否则1e-8就到头了。
3.3 sol 结构体、deval 与四输出：按需选返回形式
只要一个输出参数，ode45返回的是结构体sol，里面存了求解器内部的稠密输出信息，可以用deval在任意时刻插值求值，精度和求解器内部一致：
sol = ode45(odefun, [1 4], -2); tt = linspace(1, 4, 200); yy = deval(sol, tt); % 稠密输出插值，比 interp1 更贴求解器真实轨迹 plot(tt, yy);这一点值得多说两句。拿到[T,Y]之后用interp1做插值，本质是对已经稀疏化的采样点做二次加工，T本身就不密的地方插出来的值跟真实轨迹有偏差。deval用的是求解器在每一步内部构造的插值多项式，虽然也是近似，但它是求解器误差控制的一部分，一致性更好。做论文图、要做精细对比的时候，优先走deval。
需要事件检测时用四个输出：[T, Y, TE, YE, IE] = ode45(...)。TE是事件触发时刻，YE是对应时刻的状态，IE是触发的是第几个事件（事件函数可以返回向量，同时监控多个条件）。不设Events的时候，后三个返回空数组，写了也白写。事件函数的具体写法放到最后一章。
顺带提一句sol的另一个用处：把上次算到一半的结果当初始条件接着算。sol.x(end)和sol.y(:,end)取出末端时刻和状态，重组成新的tspan和y0再调一次，比分段手算省事。长时间积分分段跑还有个好处，中途可以检查结果是否发散，不至于跑了一夜才发现早就炸了。
4. 高阶方程降阶：把 y'' 拆成两个一阶方程再喂进去
4.1 变量替换的机械流程
ode45只认一阶方程组，二阶及以上必须先降阶。这一步没有技巧，全是固定套路，照着做就不会错。
给定 n 阶方程y⁽ⁿ⁾ = g(t, y, y', ..., y⁽ⁿ⁻¹⁾)，令x1 = y，x2 = y'，一直到xn = y⁽ⁿ⁻¹⁾。求导后得到一阶方程组：
x1' = x2 x2' = x3 ... x(n-1)' = xn xn' = g(t, x1, x2, ..., xn)初始条件也要跟着换：y0 = [y(t0); y'(t0); ...; y⁽ⁿ⁻¹⁾(t0)]。原来二阶方程给两个初值，降阶后状态向量就是二维，一一对应。
拿资料里的方程y'' = -t*y + e^t*y' + 3*sin(2t)走一遍。令x1 = y，x2 = y'，则x1' = x2，x2' = -t*x1 + e^t*x2 + 3*sin(2t)。初值y(3.9) = 2、y'(3.9) = 8，对应的y0就是[2; 8]。
用表格把对应关系固定下来，写代码时照着填：
| 原方程里的量 | 状态向量分量 | 初值 |
|---|---|---|
y | x(1) | y0(1) |
y' | x(2) | y0(2) |
y'' | 由方程组第二行给出 | 不需要 |
4.2 列向量 odefun 的写法与索引对齐
function dx = odefun(t, x) % x(1) = y, x(2) = y' dx = zeros(2, 1); % 必须预分配成列向量 dx(1) = x(2); dx(2) = -t*x(1) + exp(t)*x(2) + 3*sin(2*t); endzeros(2,1)这行是必须的。不预分配直接写dx(1) = ...，Matlab 会动态扩数组，每次扩都重新分配内存，短积分看不出来，长积分就是性能黑洞。更要命的是如果不预分配、又漏写了某个分量，返回的dx长度比y0短，ode45会报维度不匹配，报错信息指向求解器内部，看不出问题出在你的函数里。
调用：
tspan = [3.9 4.0]; y0 = [2; 8]; [t, x] = ode45(@odefun, tspan, y0); plot(t, x(:,1), '-o', t, x(:,2), '-*'); legend('y', 'y'''); xlabel('t'); ylabel('y');x(:,1)是y，x(:,2)是y'，画图时别弄反。tspan只给了 0.1 的长度，ode45会自动选步长，输出的T可能只有几个点，画出来的折线很粗，想要平滑曲线就在tspan里给一串等距时间点，或者把结果丢给deval。
阶数再高也是同一套流程。三阶方程令x1 = y、x2 = y'、x3 = y''，odefun里dx就是 3×1。新手常犯的错是“降阶降了一半”，dx(1) = x(2)写了，dx(2) = x(3)忘了，直接把y''的表达式塞进dx(2)，结果状态向量长度对不上，或者方程物理意义全错。写完odefun后默念一遍：每一行dx(i)的右边，只能出现x的各分量和t，出现y就是没换干净。
4.3 三方程耦合系统：定步长 RK4 与 ode45 的正面对照
单个方程降阶只是热身，真正体现向量化价值的是耦合方程组。资料里给的这组：
function dy = test_fun(t, y) dy = zeros(3, 1); dy(1) = y(2) * y(3); dy(2) = -y(1) + y(3); dy(3) = -0.51 * y(1) * y(2); end初值[1; 1; 3]，区间[0 15]。三个方程通过y(1)、y(2)、y(3)互相耦合，没有任何一个能单独解出来，正好用来测向量化的runge_kutta1：
[T, F] = ode45(@test_fun, [0 15], [1; 1; 3]); [T1, F1] = runge_kutta1(@test_fun, [1; 1; 3], 0.25, 0, 15); subplot(1,2,1); plot(T, F); title('ode45 自适应步长'); subplot(1,2,2); plot(T1, F1.'); title('自编 RK4, h = 0.25');注意F1.'这个转置。runge_kutta1返回的y是3×(n+1)，每列一个时刻，而plot按行画线，不转置画出来是三条竖线加一堆横线，看着像乱码。F1.'是普通转置，不是共轭转置'，这里全是实数所以两者等价，但养成用.'的习惯没坏处。
两张图放在一起，前几秒几乎重合，到t = 10之后开始分叉。原因不是 RK4 算错了，而是h = 0.25在这个系统里太粗。这组方程有正 Lyapunov 指数式的敏感依赖，初值或步长的微小差异会被指数放大，定步长的累积误差在长时间积分后压不住。要公平比较，把ode45的MaxStep也设成0.25，让它和自编版本走在同一批节点上：
opts = odeset('MaxStep', 0.25, 'RelTol', 1e-6, 'AbsTol', 1e-8); [T2, F2] = ode45(@test_fun, [0 15], [1; 1; 3], opts);这时候两者差异会小很多，剩下的差别来自ode45内部的五阶误差估计和自适应调整，属于算法本身的差异，不是实现错误。这个对照实验的价值在于：它告诉你自编函数什么时候够用、什么时候必须让位给内置求解器。短区间、平滑解、精度要求不高，自己写的定步长版本完全够；长区间、有敏感依赖、或者要求误差可控，ode45的自适应机制省下来的调试时间远超那点性能开销。
再补一个实操细节：runge_kutta1里n = floor((b-a)/h)，这里(15-0)/0.25 = 60，正好整除，不会丢尾巴。如果h取0.26，n = 57，实际只积到14.82，末端那一段凭空消失，画图时T1和T的长度不一样，plot会报错。用定步长函数之前先算一下步数，或者改成linspace生成时间轴再反推h。
5. 跑飞、报错、结果发散：五条踩坑记录与排查路径
5.1 几分钟没动静，最后抛步长警告
现象：命令窗口刷出Warning: Unable to meet integration tolerances without reducing the step size below the smallest value allowed，或者干脆卡在那里，进度条纹丝不动。
原因：方程是刚性的。显式 Runge-Kutta 的稳定域有限，遇到衰减极快的模态时，步长被稳定性而不是精度按在地上，越压越小，最后撞上浮点下限。资料里那句“若长时间没结果，应该就是刚性的”说的就是这个。
解决：换隐式求解器，ode15s是首选，ode23s、ode23tb备选。有解析 Jacobian 就通过odeset('Jacobian', @J)传进去，收敛会快不少。判断是否刚性有个土办法：先跑ode45看它是不是在某个时间段步长骤降，如果是，那段就是刚性区。
5.2 维度对不上，赋值报错
现象：Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 3-by-1，或者Index exceeds matrix dimensions。
原因：三种情况最常见。y0传了行向量，odefun返回列向量；dx没预分配，动态增长时长度跑偏；odefun里索引从 0 开始写，Matlab 是 1 基。
解决：y0统一y0(:)，odefun第一行dx = zeros(n,1)，索引检查一遍有没有x(0)。改完跑一个小tspan验证，别直接上长区间。
5.3 曲线末端炸开
现象：解在最后一段指数上升或者剧烈震荡，跟解析解完全对不上，缩小tspan又正常了。
原因：可能是方程本身在这个区间不稳定，真解就发散；可能是tspan写反了，[4 1]这种降序写法ode45会报错但有些自编函数不检查；也可能是步长越过了稳定域。
解决：先确认tspan(2) > tspan(1)。再看真解是否本来就发散，用解析解或者小步长ode15s交叉验证。如果是数值不稳定，压MaxStep，或者直接换隐式求解器。别一上来就怀疑方程写错了，先排除求解器配置。
5.4 自编 RK4 和 ode45 越比越偏
现象：两条曲线前半段重合，后半段越差越远，加大ode45的容差精度也没改善多少。
原因：比较基准不对。ode45返回的T是它自己踩的节点，自编版本返回的是等距节点，两者在完全不同的时间点上采样，直接画图看起来就是两条曲线。长时间积分还有误差放大效应，初始的微小差异被指数放大。
解决：用deval(sol, T1)在自编版本的时间点上取值再比，或者给ode45设MaxStep = h对齐步长。比较之前先明确你要验证的是“实现正确性”还是“算法精度”，前者用短区间小步长，后者才需要长区间。
5.5 改了参数结果不变
现象：改了方程里的系数、或者改了odeset里的容差，重新跑，曲线一模一样。
原因：匿名函数在创建时就把外部变量值捕获了，后面改工作区里的变量，句柄里用的还是旧值。另一种情况是odeset的参数名拼错，比如把RelTol写成Reltol（这个其实能认，不区分大小写），或者写成RelTolerance，拼错的名字被静默忽略。
解决：不要在循环里改完变量复用旧句柄，每次都重新odefun = @(t,y) ...。参数多的话用额外入参odefun(t, y, p1, p2)，较新版本的ode45支持这种写法，把参数当第四个、第五个入参传进去，比在句柄里捕获干净得多。odeset写完检查一遍拼写，不确定就把opts打印出来看。
6. 让结果站得住：事件检测与步长减半验证
事件检测是我用得最多、也最容易被忽略的功能。很多问题不需要完整积分，只想知道“状态什么时候第一次穿过零”，比如碰撞时刻、开关切换点、种群模型的灭绝时间。硬扫输出点再找符号变化，精度受采样密度限制，还容易漏掉快速穿越。Events直接把这个问题交给求解器。
事件函数签名固定，返回三个值：
function [value, isterminal, direction] = myEvent(t, y) value = y(1); % 监控第一个分量，过零即触发 isterminal = 1; % 1 = 触发时终止积分；0 = 只记录不终止 direction = 0; % 0 = 两个方向都抓；1 = 只抓上升；-1 = 只抓下降 end调用时把句柄塞进odeset，四个输出里后三个就是事件结果：
opts = odeset('Events', @myEvent, 'RelTol', 1e-8); [t, y, te, ye, ie] = ode45(@test_fun, [0 15], [1; 1; 3], opts); fprintf('首次过零时刻 %.6f，对应状态 [%.4f %.4f %.4f]\n', te(1), ye(1,:));te是触发时刻，ye是对应状态，ie是第几个事件（事件函数返回向量时可以同时监控多个条件）。isterminal设成 1 之后，积分到点就停，t和y只到那个时刻为止，后面想看也看不到，所以做完整轨迹分析时要设成 0。direction设错是常见的坑，只关心上升沿却写了 0，结果抓到的第一个点是下降沿，后面所有索引都错位。
验证数值解可信度，我固定走一遍步长减半：同一套代码，h取一组值，比较末端解的二范数差，看比值是不是接近 16。接近，说明四阶收敛成立，结果可信；远小于 16，说明已经进到舍入误差主导区，再减小步长没有意义；远大于 16 或者乱跳，说明撞上了稳定性边界或者方程本身有奇点，得回头查模型。
hs = [0.2 0.1 0.05 0.025]; err = zeros(size(hs)); ref = runge_kutta1(@test_fun, [1; 1; 3], 0.01, 0, 15); % 拿小步长当参考 for k = 1:numel(hs) [~, yk] = runge_kutta1(@test_fun, [1; 1; 3], hs(k), 0, 15); err(k) = norm(yk(:,end) - ref(:,end)); end ratio = err(1:end-1) ./ err(2:end) % 理论上应接近 16这个习惯救过我好几次。有一次曲线看着挺漂亮，收敛比却是 2 左右，一查是odefun里漏了一项，降阶之后方程根本没写全，图上的光滑只是低阶方法碰巧拟合出来的假象。从那以后我每次交付数值解结果之前，deval抽点对照和步长减半验证都强制走一遍，两个都过才敢往下用。这份资料里的自编runge_kutta1和ode45对照脚本，正好可以直接拿去做这件事，改改右端函数就能套到自己的方程上。希望帮到你。
本文还有配套的精品资源，点击获取
