当难以求得微分方程的解析解时，可以求其数值解，Matlab中求微分方程数值解的函数有七个：ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb。本文讲解ode45，其实原理都是相似的。

一、函数语法

ode45函数的形式张下面的样子

<code>[T,Y] = ode45(odefun,tspan,y0,options)

odefun: 函数句柄，用于表示ODE的右侧，例如：dy/dt=f(t,y)，可以写为f。函数句柄可以是函数，也可以是匿名函数或者内联函数tspan: 时间区间，可以是一个向量，也可是一个区间（如[t0,tf]）。y0: 初始条件，可以是一个向量。options: 可选参数，常用于控制ODE45的行为，如精度、步长等。

二、示例

我们举一个自由落体的例子来讲解ode45这个函数

{

x

˙

=

v

v

˙

=

g

\begin{cases} \dot x = v \\ \dot v = g \\ \end{cases}

{ x˙=vv˙=g​

其中

x

x

x 为位置，

v

v

v 为速度，

g

g

g 为重力加速度，设置初始值均为 0

clc;clear;

% 初始状态

state0 = [0,0];

% 时间区间

t = [0,10];

% 运算，得到T为时间点，Y为对应时间点的状态

[T,Y] = ode45(@odefunc, t, state0);

% 绘图

figure(1)

plot(T,Y(:,1));

xlabel('t(s)');

ylabel('x')

figure(2)

plot(T,Y(:,2));

xlabel('t(s)');

ylabel('v')

% 写微分方程

function dy = odefunc(t, y)

g = 9.8;

dy = zeros(2,1);

% 这里也就是求出微分方程中状态的导数

dy(1) = y(2);

dy(2) = g;

end

其中在 odefunc 函数中，

y

y

y 就可以认为是状态，而我们要给出的是下一拍状态的导数。

t

t

t 是时间范围，这里我们指定为0-10秒，state0 指的是初始状态，这里全部设定为0。T 是返回的时间点，因为ode45每一拍之间的时间间隔并不相同，他会自适应，所以需要记录下时间间隔，Y 是返回的每个时间下的状态，由于我们有两个状态，所以Y有两列。

得到位置和速度的图像

三、参数

设置参数就是设置 <code>options ，这是一个结构体，

‘RelTol’: 相对精度，默认为1e-3。‘AbsTol’: 绝对精度，默认为1e-6。‘InitialStep’: 初始步长。‘MaxStep’: 最大步长。‘Events’: 事件函数，用于检测ODE的事件点，如极值点、零点等。

如果我们想自己设置的话就是

options = odeset('RelTol',1e-4,'AbsTol',1e-8);

[T,Y] = ode45(@odefunc, t, state0);