利用ode45函数求解常微分方程的数值解

1 ode45函数的原理

该求解器有变步长和定步长两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法。

ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(Δx)^5。解决的是非刚性常微分方程。

ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode15s试试。

2 经典四阶Runge-Kutta介绍

在各种龙格－库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。

令初值问题表述如下。

y' =f(t,y)，y(t0)=y0

则，对于该问题的RK4由如下方程给出：

 y(n+1)=y(n)+h/6*(k1+2*k2+2*k3+k4)

其中

k1=f(t(n),y(n))

k2=f(t(n)+h/2,y(n)+h/2*k1)

k3=f(t(n)+h/2,y(n)+h/2*k2)

k4=f(t(n)+h,y(n)+h*k3)

这样，下一个值（yn+1）由现在的值（yn）加上时间间隔（h）和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均：

k1是时间段开始时的斜率；

k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；

k3也是中点的斜率，但是这次采用斜率k2决定y值；

k4是时间段终点的斜率，其y值用k3决定。

当四个斜率取平均时，中点的斜率有更大的权值：

RK4法是四阶方法，也就是说每步的误差是h阶，而总积累误差为h阶。

3 ode45函数求解常微分方程举例

求解微分方程d²y/dt²- 7(1-y²)dy/dt + y = 0在初始条件下y(0) = 1、y’(0) = 0下的解，并画出解的图。

解：设x1 = y，x2 = dy/dt，将二阶方程化成一阶方程组

   dx1 = x2  x1(0) = 1;

   dx1/dt = 7(1-x1²)*x2-x1  x2(0) = 0;

在北太天元脚本编辑器窗口输入以下程序，并保存文件名为vdp.m。

function fy = vdp(t,x)

fy = [x(2);7*(1-x(1)^2)*x(2)-x(1)];

在北太天元命令行窗口输入以下程序。

>> y0 = [1;0];

>> [t,x] = ode45(@vdp,[0 40],y0);

>> y = x(:,1);

>> plot(t,y)

>> xlabel('t')

>> ylabel('y')

输出结果如下图所示：