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函數求解常微分方程舉例

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

解：設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')

輸出結果如下圖所示：