简介：围绕MATLAB求解常微分方程初值问题的核心函数ode45，这份资料系统性整理了函数调用格式、dydt方程定义、参数设置、指定输出点、事件检测、多输出系统等关键用法，并配有可运行的.m示例脚本。ode45基于经典四阶龙格-库塔方法，适用于非线性、高阶及随时间变化的微分方程组，适合需要掌握该求解器的本科生、研究生和工程技术人员。压缩包共4个文件，包含2个PDF讲义、1个PPT课件和1个m脚本，总大小14.65MB；两份PDF分别讲解ode函数与其他函数的配合使用以及MATLAB/Simulink中的ode45程序实现，PPT则结合动力学与振动工程案例展示实际建模与求解过程。已有2017人学习下载，内容兼顾原理梳理与动手实践，既有细致笔记和可执行脚本，也有工程课件，可作为初学者进阶和日常仿真的常备参考。
1. 为什么 MATLAB 求解微分方程总绕不开 ode45
做动力学与振动仿真的人，多半在《MATLAB与工程应用》第 7 章讲义里第一次碰到 ode45：给定 dy/dt=f(t,y) 与初值 y(t0)，求出状态随时间的变化。这份《ode45的笔记》收录了 ode45xuexi.m 示例脚本、ode 使用数据库 PDF、动力学与振动 PPT，以及 MATLAB 与 Simulink 双版本用法文档，覆盖从调用到后处理的完整链路。
对刚接触 matlab 求解微分方程的人，ode45 是第一个该掌握的求解器，默认参数就能跑通大部分非线性、高阶、时变系统；对写过几年仿真代码的人来说，真正值得抠的是容差、参数透传、事件检测这些边角，笔记里都有对应内容。
下文按原理、调用、参数、高阶与刚性、验证展开，代码可直接运行，建议开着 MATLAB 边看边敲。
2. Runge-Kutta 45 原理与 ode45 函数签名
2.1 Dormand–Prince 自适应步长：ode45 的“45”指什么
ode45 底层是 Dormand–Prince 对（RK45），本质是一对显式 Runge-Kutta 公式：同一个积分步里同时算出一个四阶近似和一个五阶近似，两者之差就是局部截断误差的估计，求解器据此决定下一步长放大还是缩小。这个自适应机制让 ode45 每一步尽量取大步长，又能把误差压在容差内，这正是它被 MATLAB 官方列为默认首选求解器的原因。
这个设计带来两个使用结论。一是 ode45 名义上是四阶方法，但步长控制用的是五阶值，实际精度介于四阶与五阶之间；二是它只适合非刚性（non-stiff）问题，显式公式的稳定域有限，遇到时间常数相差几个数量级的系统，步长会被压得非常小，届时需要切到隐式求解器。
提示：官方文档 Choosing a Solver 把 ode45 列为大多数问题的 first try，但它的定位是“均衡”而不是“高精度”。容差要求极严或方程刚硬时，要换 ode113、ode15s。
2.2 最小可用调用：tspan、y0 与函数句柄
先跑通最简形式。以指数衰减方程 dy/dt=-0.5y、初值 y(0)=1 为例：
f = @(t, y) -0.5 .* y; % 匿名函数定义方程右侧 [t, y] = ode45(f, [0 10], 1); plot(t, y, 'o-'); xlabel('t'); ylabel('y'); grid on;f 是函数句柄，MATLAB 规定它必须先收 t 再收 y，即使方程与 t 无关也要保留第一个形参，返回值就是 dy/dt。tspan 给两个元素 [0 10] 表示只声明积分区间，不关心中间输出点，步长完全交给自适应机制；y0 是 tspan(1) 处的初始状态，标量方程给标量，多变量系统给列向量。返回的 t 是算法自动选出的时间点，y 是与 t 等长的解序列，曲线很密的地方说明步长被容差逼小了，很疏的地方暗示步长已经顶到上限。
ode45 的完整签名如下，工程代码里多数情况只用前四个参数：
| 参数 | 含义 | 典型写法 |
|---|---|---|
| odefun | 方程右侧函数句柄，接收 (t,y) 返回 dy/dt | @(t,y) -0.5*y、@myode |
| tspan | 积分区间，或带输出点要求的向量 | [0 10]、0:0.01:10 |
| y0 | 初始状态，维度与 dy/dt 返回量一致 | 1、[0; 0] |
| options | odeset 生成的选项结构体 | odeset('RelTol',1e-6) |
签名里 options 之后其实还留有第五、第六个附加参数位，老式写法会往里塞透传变量，比如ode45(@dydt,tspan,y0,[],[],k)：空方括号分别是 options 和第一个透传参数的占位，k 透传给 dy/dt 作为附加输入，因此 dydt 的形参必须按 (t,y,p1,p2) 对齐，个数对不上照样报错。这个写法在兼容模式下能跑，但可读性和健壮性都不如用匿名函数或嵌套函数绑定参数。
2.3 在 matlab 中定义微分方程：三种写法与参数绑定
工程模型很少只有一个方程，这里用经典的 van der Pol 方程展示三种定义方式，它本身常被写成二阶系统降阶后的两变量方程组：
% 方法一：独立函数文件，适合方程长、规模大的模型 function dydt = vdp_ode(t, y) mu = 1; % 参数写在函数体里，固定 dydt = [y(2); mu*(1-y(1)^2)*y(2) - y(1)]; end % 主程序调用： [t, y] = ode45(@vdp_ode, [0 20], [2; 0]); % 方法二：匿名函数绑定工作区变量，参数扫描最常用 mu = 2; f = @(t, y) [y(2); mu*(1-y(1)^2)*y(2) - y(1)]; [t, y] = ode45(f, [0 20], [2; 0]); % 方法三：嵌套函数共享主函数局部变量，脚本更干净 function run_vdp() mu = 2; % 嵌套函数能直接看到 mu [t, y] = ode45(@dydt, [0 20], [2; 0]); plot(t, y(:,1)); function dydt = dydt(t, y) dydt = [y(2); mu*(1-y(1)^2)*y(2) - y(1)]; end end三种写法对应三种场景。方法一是模型文件里最常见的组织方式，参数固定、接口清晰，多个算例复用时不用复制方程体；方法二把 mu 留在当前工作区，改参数后重新构造句柄即可，做参数扫描循环时最省事，但循环里要先把当前 mu 复制到局部变量再构造句柄，否则句柄捕获的是循环退出后的最终值；方法三把辅助函数藏进主函数内部，脚本文件结构完整，适合一个脚本解决一个完整算例。y 的两个分量 y(1)、y(2) 对应状态变量，dydt 必须返回列向量，行向量写法在某些版本会静默转置并给警告，建议始终显式写成分号分隔。
3. odeset 容差参数与输出点控制
3.1 RelTol 与 AbsTol：误差公式与取值策略
ode45 的步长控制目标是让每一步的局部误差 e(i) 满足不等式 |e(i)| ≤ max(RelTol·|y(i)|, AbsTol(i))。RelTol 按当前状态量级缩放误差上限，AbsTol 给绝对底线，二者取较大值。这样既不会因为状态量级很大而要求离谱的绝对精度，也不会因为状态接近 0 而让相对误差失去意义。
默认 RelTol=1e-3、AbsTol=1e-6，画曲线形态够用；如果要拿结果做微分、积分、参数辨识，建议收敛到 1e-6 至 1e-9。当某个状态分量长期贴着 0 走，比如振动分析里的平衡位置附近，相对误差会失效，必须把对应分量的 AbsTol 收紧，否则步长控制会失灵：
opts = odeset('RelTol', 1e-6, 'AbsTol', [1e-8 1e-8]); [t, y] = ode45(@dydt, [0 30], [0; 0], opts);options 是第四个位置参数，AbsTol 给向量时长度与状态数一致，给标量时所有分量共用。odeset 支持一次设置多个字段，未设置的项沿用默认。ode 使用与其他函数使用数据库 PDF 里整理了完整参数表，实际项目里高频用到的是这几个：
| 选项 | 默认值 | 作用 |
|---|---|---|
| RelTol | 1e-3 | 相对误差容限，决定解的相对精度 |
| AbsTol | 1e-6 | 绝对误差容限，决定接近零分量的精度 |
| MaxStep | 区间长度的 1/10 | 限制最大步长，防跳过事件或控制计算量 |
| InitialStep | 算法估算 | 手动给定第一步长，用于强瞬态问题 |
| Stats | off | 置 'on' 打印成功与失败步数 |
| Events | 无 | 事件函数句柄，见第 4 章 |
3.2 固定输出点与 deval 任意取值
很多时候图纸需要等间距时间轴的数据，两种做法效果差别很大：
% 方式一：tspan 给稠密向量，输出点严格落在指定坐标 tspan = 0:0.05:30; [t, y] = ode45(@dydt, tspan, [0; 0]); % 方式二：一次自适应积分，之后任意位置插值 sol = ode45(@dydt, [0 30], [0; 0]); tq = linspace(0, 30, 500); yq = deval(sol, tq);方式一里列出的每个点都会进入求解器内部的自适应判定流程，点数设得太密会拖慢求解。方式二只在 [0 30] 上做一次积分，deval 在每步构造的多项式插值上取任意坐标，适合积分一次、多处取数。这份笔记里出现的ode45(@dydt,tspan,y0,[],[],tout)是早期教程的写法，实际效果是把 tout 作为附加参数透传给 odefun，对输出点并没有强制约束力；要固定输出点，正路是上面两种。绘图和后续数组运算用方式一更直白，需要高分辨率曲线或对局部区间加密时，方式二更省计算。
3.3 用 Stats 看求解器是否在硬扛
性能问题要先看数字而不是猜：
opts = odeset('Stats', 'on'); [t, y] = ode45(@dydt, [0 30], [0; 0], opts);运行后命令行输出 successful steps 与 failed attempts 的计数。failed attempts 占比明显偏高，说明步长频繁被拒，要么容差给得太严，要么系统已经接近刚性。把这个数字和 MaxStep 的设置放在一起看，是定位“ode45 为什么这么慢”的第一步。
4. 多变量动力学系统、事件检测与刚性方程切换
4.1 二阶系统降阶：质量-弹簧-阻尼的状态空间写法
《MATLAB与工程应用》第 7 章的动力学与振动问题里，最典型的是 mx''+cx'+kx=F(t)，这也是那份 PPT 里反复出现的模型。二阶方程必须先降阶才能交给 ode45：令 y1=x、y2=x'，得到 dy1/dt=y2、dy2/dt=(F(t)-c·y2-k·y1)/m。在 matlab 中定义微分方程时，多变量系统就是把所有导数写成一个列向量：
m = 2; c = 0.6; k = 8; F0 = 1.5; w = 2; % 正弦激励幅值与角频率 F = @(t) F0 * sin(w * t); % 激励也可换成阶跃、脉冲或任意自定义函数 dydt = @(t, y) [y(2); (F(t) - c*y(2) - k*y(1)) / m]; [t, y] = ode45(dydt, [0 30], [0; 0]); plot(t, y(:,1), t, y(:,2)); legend('位移 x', '速度 v'); grid on;dydt 返回两行列向量，第一行是位移的导数即速度，第二行是加速度的显式表达式；y0=[0;0] 表示零初始位移、零初速，长度必须与 dydt 返回的列数一致。输出的 y 是 N×2 矩阵，第 1 列位移、第 2 列速度，后续算机械能、画相轨迹、做频谱分析都从这两列取数。外力 F(t) 写成句柄后可以随时替换为阶跃函数、矩形脉冲或实测数据插值，这是 ode45 处理时变系统最方便的地方。
4.2 事件检测：落地时刻与过零触发
ode45 按时间积分，本身不知道“什么时候撞到地面”。Events 机制监控某个量的过零，一旦符号变化就记录时刻，并可选终止积分。以自由落体问落地时间为例：
function dy = fall_ode(t, y) g = 9.8; dy = [y(2); -g]; % y(1) 高度，y(2) 下落速度 end function [value, isterminal, direction] = hit_ground(t, y) value = y(1); % 监控高度 isterminal = 1; % 触发后停止积分 direction = -1; % 只在高度由正变负时触发 end opts = odeset('Events', @hit_ground, 'RelTol', 1e-8); [t, y, te, ye] = ode45(@fall_ode, [0 10], [10; 0], opts); fprintf('落地时刻 t=%.6f，速度=%.6f\n', te(end), ye(end,2));事件函数必须返回三个量。value 是被监控量，其符号变化触发事件；isterminal 取 1 则触发后终止，取 0 则记录每次过零继续积分，适合多次碰撞或周期事件；direction 取 1、-1、0 分别限定正向过零、负向过零、双向触发。注意如果 value 在积分区间内一次都没过零，te 返回空数组，直接索引 end 会报错，先判断 isempty(te)。
4.3 刚性判定与 ode15s 切换
刚性的直观表现是系统里同时存在快、慢两类时间常数。显式 RK45 为保证稳定性，步长必须小于最快分量的尺度，于是总步数爆炸。三个诊断信号：Stats 的 failed attempts 异常偏高、MaxStep 被算法自动压到极小、相同区间 ode45 耗时远超 ode23 或 ode15s。用一个经典刚性方程做实验：
f_stiff = @(t, y) -1000*(y - sin(t)) + cos(t); % 快衰减叠加慢正弦 tic; [t1, y1] = ode45(f_stiff, [0 10], 1); toc; tic; [t2, y2] = ode15s(f_stiff, [0 10], 1); toc; fprintf('ode45 步数 %d，ode15s 步数 %d\n', numel(t1), numel(t2));老版本 MATLAB 上这个对比非常悬殊，ode45 上万步、ode15s 只有几十步；如果你的版本做了隐式刚度检测优化，对比不明显，把系数 1000 改成 1e5 再观察耗时差距。求解器选型按下表：
| 求解器 | 方法 | 适用场景 |
|---|---|---|
| ode45 | 显式 RK45 | 大多数非刚性问题首选 |
| ode23 | 显式 RK23 | 容差要求不高、积分区间较短 |
| ode113 | 变阶 Adams | 非刚性、容差很严、被积函数平滑 |
| ode15s | 变阶隐式 | 刚性问题，或 ode45 计算量不可接受 |
| ode23s | 低阶隐式 | 刚性但容差要求粗、希望步长尽量大 |
Simulink 模型的默认求解器同样是 ode45，当模型里存在惯性相差很大的环节，比如快变电气量与慢变机械量耦合，要在 Solver 配置面板里手动切 ode15s，否则仿真时长会指数级上涨，这是《ode45的用法和程序（matlab和simulink）》那份 PDF 里专门提醒过的一处。
5. 收敛性验证与 ode45 步长陷阱
5.1 解析解对拍
笔记里的 ode45xuexi.m 把调用、事件、验证串了一遍，演示脚本最后留下来的检查习惯是这三条。第一条，拿到数值解先不急着画图，用已知解析解做一次收敛性检查。以 dy/dt=-y、y(0)=1 为例，解析解是 exp(-t)：
t_ana = linspace(0, 5, 100); y_ana = exp(-t_ana); for rt = [1e-3, 1e-6, 1e-9] opts = odeset('RelTol', rt); [t, y] = ode45(@(t, y) -y, [0 5], 1, opts); err = max(abs(interp1(t, y, t_ana) - y_ana)); fprintf('RelTol=%g, 最大误差=%g\n', rt, err); endRelTol 从 1e-3 收紧到 1e-9，最大误差应单调下降两到三个量级；如果误差不降反升或出现平台，先检查方程有没有写反，而不是怀疑求解器。interp1 的作用是把不同步长的数值解对齐到同一组采样坐标，比较才有意义。
5.2 用不变量检查长期积分漂移
第二条，对振荡系统盯不变量。无阻尼弹簧振子的机械能 E=0.5·m·y(2)²+0.5·k·y(1)² 应为常数，积分误差累积会让它缓慢漂移：
E = 0.5*m*y(:,2).^2 + 0.5*k*y(:,1).^2; drift = abs(E(1) - E(end)) / E(1); if drift > 1e-6 opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10); [t, y] = ode45(dydt, [0 30], [0; 0], opts); end能量漂移是长期积分可信度的直接度量。振荡问题局部误差虽小，长时间积分后趋势性漂移会被放大，只看波形重合不够，不变量检查能暴露系统性偏差。
5.3 四个高频翻车点
第三条，把常见误用记成清单。tspan 写成 [10 0] 倒序，MATLAB 仍会积分，但输出时间轴反向，绘图与后续插值全部错位，统一用升序区间；dydt 返回行向量 [a, b] 而不是 [a; b]，部分版本会静默转置并给警告，换版本行为可能不一致，一律显式返回列向量；y0 维度与 dydt 返回长度不一致时，报错指向数组索引越界，先核对状态变量个数再查 tspan；MaxStep 设置不当，过大可能直接跨过事件触发点导致 te 为空，过小则让输出数组膨胀到内存吃紧。我现在的习惯是先开 Stats 看默认步数，再按步数的十分之一量级设定 MaxStep，兼顾事件捕获与资源开销。
本文还有配套的精品资源，点击获取
