matlab数值计算函数--ode45

点击全文阅读

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

一、函数语法

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

[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有两列。

得到位置和速度的图像

三、参数

设置参数就是设置 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);

点击全文阅读