在使用 ODE 指令時,我們必須先將要求解的 ODE 表示成一個函式,此函式的輸入為 t(時間)及 y(狀態變數,State Variables),輸出則為 dy(狀態變數的微分值)。如果此 ODE 函式(或又稱 ODE 檔案)的檔名為 odeFile.m,則我們呼叫 ODE 指令的格式如下:
[t, y] = solver ('odeFile', [t0, t1], y0); 其中 [t0, t1] 是積分的時間區間,y0 代表起始條件(Initial Conditions),solver 是前表所列的各種 ODE 指令,t 是輸出的時間向量,y 則是對應的狀態變數向量。
以 van der Pol 微分方程為例,其方程式為:
$$y''-\mu (1-y^2)y'+y=0$$首先我們必須將它化成標準格式,令 $y_1=y$,$y_2=y'$,則上述微分方程式可以改寫為兩個一階的微分方程式:
$$ \left\{ \begin{array}{rcl} y_1' & = & y_2 \\ y_2' & = & \mu (1-y_1^2)y_2-y_1\\ \end{array} \right. $$此種標準格式可用向量來表示成一般化的數學式:
$$\vec{y}' = \vec{F}(\vec{y}, t)$$其中 $\vec{y}'$ 為一向量,代表狀態變數。使用此標準格式後,並假設 $\mu=1$,我們的 ODE 檔案(vdp1.m)可顯示如下:
一旦有了 ODE 檔案,我們即可選用前述之 ODE 指令來朮解。在 $\mu=1$ 時,van der Pol 方程式並非 Stiff 系統,因此,我們可選用 ode45 來求解,最簡單的作法,就是使用此指令來畫出積分的結果,範例如下:
其中 [0, 25] 代表積分的時間區間,[3 3]’ 則代表起始條件(必須以行向量來表示)。由於沒有輸出變數,所以上述程式執行結束後,MATLAB 只會畫出狀態變數對時間的圖形,如上圖。若要取得積分所得的狀態變數及對應的時間,可以加上輸出變數以取得這些資料,範例如下:
我們也可以畫出 $y(t)$ 及 $y'(t)$ 在相位平面(Phase Plane )的運動情況,例如:
在上例中,當 $\mu$ 值越來越大時,van der Pol 方程式就漸漸變成一個 Stiff 的系統,此時若要解此方程式,就必須改用專門對付 Stiff 系統的指令。
首先,我們將 $\mu$ 值改成 1000, ODE 檔案要改寫如下(vdp2.m):
此時即可選用專門對付 Stiff 系統的 ODE 指令,例如 ode15s,來求解此系統並作圖顯示,如下:
由上圖可看出,$y'(t)$ 的變化相當劇烈(超過 $\pm 1000$),這就是 Stiff 系統的特色。
若要畫出二維平面相位圖,可以使用下列範例:
若要產生在某些特定時間點的狀態變數值,則我們呼叫 ODE 指令的格式可改成:
[t, y] = solver('odeFile', [t0, t1, ..., tn], y0); 其中 [t0, t1, ..., tn] 即是特定時間點所形成的向量,請讀者自行試試看。
MATLAB程式設計:進階篇