1 ode45函數(shù)的原理

該求解器有變步長(zhǎng)和定步長(zhǎng)兩種類型。不同類型有著不同的求解器，其中ode45求解器屬于變步長(zhǎng)的一種，采用Runge-Kutta算法。

ode45表示采用四階-五階Runge-Kutta算法，它用4階方法提供候選解，5階方法控制誤差，是一種自適應(yīng)步長(zhǎng)（變步長(zhǎng)）的常微分方程數(shù)值解法，其整體截?cái)嗾`差為(Δx)^5。解決的是非剛性常微分方程。

ode45是解決數(shù)值解問(wèn)題的首選方法，若長(zhǎng)時(shí)間沒(méi)結(jié)果，應(yīng)該就是剛性的，可換用ode15s試試。

2 經(jīng)典四階Runge-Kutta介紹

在各種龍格－庫(kù)塔法當(dāng)中有一個(gè)方法十分常用，以至于經(jīng)常被稱為“RK4”或者就是“龍格－庫(kù)塔法”。該方法主要是在已知方程導(dǎo)數(shù)和初值信息，利用計(jì)算機(jī)仿真時(shí)應(yīng)用，省去求解微分方程的復(fù)雜過(guò)程。

令初值問(wèn)題表述如下。

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

則，對(duì)于該問(wèn)題的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)

這樣，下一個(gè)值（yn+1）由現(xiàn)在的值（yn）加上時(shí)間間隔（h）和一個(gè)估算的斜率的乘積所決定。該斜率是以下斜率的加權(quán)平均：

k1是時(shí)間段開(kāi)始時(shí)的斜率；

k2是時(shí)間段中點(diǎn)的斜率，通過(guò)歐拉法采用斜率k1來(lái)決定y在點(diǎn)tn+h/2的值；

k3也是中點(diǎn)的斜率，但是這次采用斜率k2決定y值；

k4是時(shí)間段終點(diǎn)的斜率，其y值用k3決定。

當(dāng)四個(gè)斜率取平均時(shí)，中點(diǎn)的斜率有更大的權(quán)值：

RK4法是四階方法，也就是說(shuō)每步的誤差是h階，而總積累誤差為h階。

3 ode45函數(shù)求解常微分方程舉例

求解微分方程d2y/dt2- 7(1-y2)dy/dt + y = 0在初始條件下y(0) = 1、y’(0) = 0下的解，并畫(huà)出解的圖。

解：設(shè)x1 = y，x2 = dy/dt，將二階方程化成一階方程組

   dx1 = x2  x1(0) = 1;

   dx1/dt = 7(1-x12)*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')

輸出結(jié)果如下圖所示：