11-2 ODE O禰誑峈k

在使用 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)可顯示如下:

11-常微分方程式/vdp1.mfunction dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

一旦有了 ODE 檔案,我們即可選用前述之 ODE 指令來朮解。在 $\mu=1$ 時,van der Pol 方程式並非 Stiff 系統,因此,我們可選用 ode45 來求解,最簡單的作法,就是使用此指令來畫出積分的結果,範例如下:

Example 1: 11-常微分方程式/odeBasic01.mode45('vdp1', [0 25], [3 3]');

其中 [0, 25] 代表積分的時間區間,[3 3]’ 則代表起始條件(必須以行向量來表示)。由於沒有輸出變數,所以上述程式執行結束後,MATLAB 只會畫出狀態變數對時間的圖形,如上圖。若要取得積分所得的狀態變數及對應的時間,可以加上輸出變數以取得這些資料,範例如下:

Example 2: 11-常微分方程式/odeGetData01.m[t, y] = ode45('vdp1', [0 25], [3 3]'); plot(t, y(:,1), t, y(:,2)); xlabel('Time t'); ylabel('Solution y(t) and y''(t)'); legend('y(t)', 'y''(t)');

我們也可以畫出 $y(t)$ 及 $y'(t)$ 在相位平面(Phase Plane )的運動情況,例如:

Example 3: 11-常微分方程式/odephasePlot01.m[t, y] = ode45('vdp1', [0 25], [3 3]'); plot(y(:,1), y(:,2), '-o'); xlabel('y(t)'); ylabel('y''(t)');

在上例中,當 $\mu$ 值越來越大時,van der Pol 方程式就漸漸變成一個 Stiff 的系統,此時若要解此方程式,就必須改用專門對付 Stiff 系統的指令。

首先,我們將 $\mu$ 值改成 1000, ODE 檔案要改寫如下(vdp2.m):

11-常微分方程式/vdp2.mfunction dy = vdp2(t, y) mu = 1000; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

此時即可選用專門對付 Stiff 系統的 ODE 指令,例如 ode15s,來求解此系統並作圖顯示,如下:

Example 4: 11-常微分方程式/ode15s01.m[t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(2,1,1); plot(t, y(:,1), '-o'); xlabel('Time t'); ylabel('y(t)'); subplot(2,1,2); plot(t, y(:,2), '-o'); xlabel('Time t'); ylabel('y''(t)'); % 注意單引號「'」的重覆使用

由上圖可看出,$y'(t)$ 的變化相當劇烈(超過 $\pm 1000$),這就是 Stiff 系統的特色。

若要畫出二維平面相位圖,可以使用下列範例:

Example 5: 11-常微分方程式/ode15s02.m[t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(1,1,1); plot(y(:, 1), y(:, 2), '-o'); xlabel('y(t)'); ylabel('y''(t)')

若要產生在某些特定時間點的狀態變數值,則我們呼叫 ODE 指令的格式可改成:

[t, y] = solver('odeFile', [t0, t1, ..., tn], y0);

其中 [t0, t1, ..., tn] 即是特定時間點所形成的向量,請讀者自行試試看。


MATLAB程式設計:進階篇