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

标签: 科学计算

搬砖的攻城狮 2023-03-09 09:17:17

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')

 

输出结果如下图所示:

 

微信截图_20230309091828.png

1980 0 1 收藏 回复

回复

回复

重置 提交