matlab数值计算函数--ode45

当难以求得微分方程的解析解时，index.php/tags-653.html" class="superseo">可以求其数值解，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);