第二十九章

第二十九章

常微分方程式

本 章 重 點

本 章 介 紹 MATLAB 用 於 求 解 常 微 分 方 程 式

（Ordinary Differential Equations，簡稱 ODE）的指

令，這些指令可對一般的 ODE 或 Stiff 的 ODE 進

行數值積分求解。

29-1 ODE 指令簡介

MATLAB 用於求解常微分方程式的指令可以列表如下：

整理：

指令 方法 說明

ode45 Explicit Runge-Kutta（4, 5）pair of Dormand-Prince Nonstiff Solver

ode23

Explicit Runge-Kutta （2, 3）pair of Bogacki and Shampine

Nonstiff Solver

ode113 Variable order Adams-Bashforth-Moulton PECE solver

Nonstiff Solver ode15s Numerical differentiation formulas（NDFS） Stiff Solver ode23s Modified Rosenbrock formula of order 2 Stiff Solver ode23t Trapezoidal rule with a“free”interpolant Stiff Solver

ode23tb

Implicit Runge-Kutta formula with a backward

differentiation formula of order two Stiff Solver

提 示 ：

8 除了 ODE 指令外，您也可使用 Simulink 來求解常微分方程式。

上表列出的指令項目繁多，但最主要可分兩大類，即適用於 Nonstiff 系 統及 Stiff 系統的指令。一般的常微分方程式都是 Nonstiff 系統，所以 可以直接選用 ode45、ode23 或 ode113 來求解。Stiff 系統的特性是它 的速率（即微分值）差異相常大，因此若使用一般的 ode45、ode23 或

ode113 來求解，可能會使得積分的步長（Step Sizes）變得很小（以降 低積分誤差至可容忍範圍以內），導致計算時間過長。因此，欲解決 此類 Stiff 系統 問題， 只要選 用專門 對付 Stiff 系統 的指令 ，例 如 ode15s、ode23s、ode23t 及 ode 23tb，即可迎刃而解。

29-2 ODE 指令基本用法

在使用 ODE 指令時，我們必須先將要求解的 ODE 表示成一個函式，

此函式的輸入為 t（時間）及 y（狀態變數），輸出則為 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 微分方程為例，其方程式為：

0 ' ) 1 (

'' − − y

y + y = y µ

首先我們必須將它化成標準格式，令

y

= y

y

= y '

式可以重寫為：

 



−

−

=

=

1 2 2 1 2

2 1

) 1 ( '

'

y y y y

y y

µ

此種標準格式可用數學式表示成：

) , ( ' F y t y r r r =

其中

y r

為一向量，代表狀態變數（State Variables）。使用此標準格式 後，並假設

µ

>> type vdp1.m

function dy = vdp1(t, y) mu = 1;

dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

一旦有了 ODE 檔案，我們即可選用 ODE 指令來朮解。在

µ

van der Pol 方程式並非 Stiff 系統，因此，我們可選用 ode45 來求解，

最簡單的作法如下：

>> ode45('vdp1', [0 25], [3 3]');

0 5 10 15 20 25

-3 -2 -1 0 1 2 3 4

其中 [0, 25] 代表積分的時間區間，[3 3] 則代表起始條件（必須以行向 量來表示）。程式執行結束後，MATLAB 只會畫出狀態變數對時間的 圖形，如上圖。若要取得積分所得的狀態變數及對應的時間，可加上 輸出變數，例如：

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

0 5 10 15 20 25

-3 -2 -1 0 1 2 3 4

Tim e t

Solution y(t) and y'(t)

y(t) y'(t)

我們也可以畫出

y (t )

y ' t ( )

>> plot(y(:,1), y(:,2), '-o');

>> xlabel('y(t)'); ylabel('y''(t)');

-3 -2 -1 0 1 2 3 4 -3

-2 -1 0 1 2 3

y(t)

y'(t)

在上例中，當

µ

首先，我們將

µ

>> type vdp2.m

function dy = vdp2(t, y) mu = 1000;

dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

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

>> [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)'); % 注意單引號「'」的重覆使用

0 500 1000 1500 2000 2500 3000 -4

-2 0 2 4

Tim e t

y(t)

0 500 1000 1500 2000 2500 3000

-2000 -1000 0 1000 2000

由上圖可知，

y ' ( t )

± 1000

>> subplot(1,1,1);

>> plot(y(:, 1), y(:, 2), '-o');

>> xlabel('y(t)'); ylabel('y''(t)');

y(t)

y'(t)

-3 -2 -1 0 1 2 3 -1500

-1000 -500 0 500 1000 1500

y(t)

y'(t)

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

[t, y] = solver('odefile', [t0, t1, … , tn], y0);

其中 [t0, t1, … , tn] 即是特定時間點所形成的向量。

29-3 ODE 指令的選項

所有的 ODE 指令可以接受第四個輸入變數，代表積分過程用到的各 種選項（Options），此種 ODE 指令的格式為：

[t, y] = solver('odefile', [t0, tn], y0, options);

其中 options 是由 odeset 指令來控制的結構變數，此結構變數即包含了 積分過程用到的各種選項，odeset 的一般格式如下：

二維相位平面圖

options = odeset('name1', value1, 'name2', value2, … )

此時產生的結構變數 options，其欄位 name1 的值為 value1，欄位 name2 的值為 value2，依此類推。未被設定的欄位，其欄位值即為預設值。

我們也可以只改變一個現存的 options 結構變數中，某個欄位的值，其 格式如下：

new_options = odeset(options, 'name', value);

若要讀出某個欄位的值，可用 odeget，其格式如下：

value = odeget(otpions, 'name');

其中 name 為欄位名稱，value 即為對應的欄位值。

當使用 odeset 指令時，若無任何輸入變數，則 odeset 會顯示所有的欄 位名稱及欄位值，並以大括號代表預設值，例如：

>> odeset

AbsTol: [ positive scalar or vector {1e-6} ] BDF: [ on | {off} ]

Events: [ on | {off} ] InitialStep: [ positive scalar ] Jacobian: [ on | {off} ] JConstant: [ on | {off} ] JPattern: [ on | {off} ]

Mass: [ {none} | M | M(t) | M(t,y) ] MassSingular: [ yes | no | {maybe} ]

MaxOrder: [ 1 | 2 | 3 | 4 | {5} ] MaxStep: [ positive scalar ] NormControl: [ on | {off} ] OutputFcn: [ string ]

OutputSel: [ vector of integers ] Refine: [ positive integer ]

RelTol: [ positive scalar {1e-3} ]

Stats: [ on | {off} ] Vectorized: [ on | {off} ]

這些欄位名稱及相關說明，可見下表：

整理：由 odeset 產生的 ODE 選項

類別 欄位名稱 資料型態 預設值 說明

RelTol 正純量

10

誤差容忍度之

相關欄位 AbsTo1 正純量或向量

10

OutPutFcn 字串 ‘odeplot’

輸出函式（若 ODE 指令 無輸出變數，則在數值積 分執行完畢後，MATLAB 會呼叫此輸出函式）

OutputSel 索引向量 全部

ODE 指令之輸出變數的 索引值，以決定那些輸出 變數之元素將被送到輸 出函式

Refine 正整數 1 或

4 (for ode45)

Refine = 2 可產生兩倍數 量的輸出點，Refine = 3 可產生三倍數量的輸出 點，依此類推。

積分輸出之相 關欄性

Stats on 或 off off

Stats = 'on' 產生計算過 程的各種統計資料。

Jacobian 矩陣

之相關欄位 Jconstant on 或 off off

如果 Jacobian 矩陣常 數，則 JConstant = 'on'

Jacobian on 或 off off

若 F(t, y, ’Jacobian') 傳

回 y F

∂

∂

，則 Jacobian = 'on'

Jpattern on 或 off off

若 F([ ], [ ], ’JPattern’) 傳

回 y F

∂

∂

，且 y F

∂

∂

是稀疏矩 陣，則 JPattem = 'on'

Vectorized on 或 off off

若 F(t, [y1, y2… ..]) 傳回 [F(t,y1), F(t,y2)… ..]，則 Vectorized = 'on'

Max Step 正純量 ODE 指令之積分步長 的

上限 積分步（Step

Sizes）之相關

欄位 Initial Step 正純量 起始步長的建議值

Mass

none，M，

M(t)，或 M(t, y) none

表明 ODE 指令案是否會 傳回質量矩陣

質量矩陣之相

關欄位 MassSingu lar

yes，no 或

maybe maybe

表明質量矩陣是否為 Singular

事件發生時間

之相關欄位 Events on 或 off off

若 ODE 檔案並傳回事件 函式及相關資訊，則 Events = 'on' Ode15s 之相

關欄位 MaxOrder 1, 2, 3, 4 或 5 5

積分公式用到的最高次 數

BDF on 或 off off

若使用 BDF（Backward Differentiation

Formula）則 BDF = 'on'

以下對幾個常用到的欄位來進行說明。

f 在積分誤差容忍度方面，每一次積分所產生的局部誤差 e(i)，必須滿 足下列方程式：

≤ (i )

e

y (i )

其中 i 代表第 i 個狀態變數，因此，我們可降低 RelTol 及 AbsTol 來求 得更精確的積分結果。例如：

>> options = odeset('RelTol', 0.0001);

>> ode45('vdp1', [0 25], [3 3]', options);

0 5 10 15 20 25

-3 -2 -1 0 1 2 3 4

在上圖中，由於相對誤差值很小，因此我們得到較細密的點，讀者可 將之與本章第二節的第一個圖比較，即可驗證。

在積分輸出的相關處理方面，我們可以選用一個 OutputFcn，使得當 ODE 指令沒有輸出變數時，此輸出函式 OutputFcn 會被 MATLAB 呼 叫。OutputFcn 的預設值是“odeplot”，其功能為畫出所有的狀態變 數，其它可用的函式還有：

l odephas2：畫出 2-D 的相位平面（Phase Plane）

l odephas3：畫出 3-D 的相位平面

l odeprint：隨時在指令視窗印出計算結果

以 Lorenz 渾沌方程式（Lorenz Chaotic Equation）為例，其 ODE 檔案 之內容可顯示如下：

>> type lorenz1.m

function dy = lorenz1(t, y)

% LORENZ1: ODE file for Lorenz chaotic equation SIGMA = 10.;

RHO = 28.;

BETA = 8./3.;

A = [ -BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1 ];

dy = A*y;

以 ode45 進行數值積分之結果可繪圖如下：

>> ode45('lorenz1', [0 10], [20 5 -5]');

0 2 4 6 8 10

-30 -20 -10 0 10 20 30 40 50

三個狀態變數

若要畫三度空間之相位圖形，可輸入如下：

>> options = odeset('OutputFcn', 'odephas3'); % 使用 odephas3

% 進行繪圖

>> ode45('lorenz1', [0 10], [20 5 -5]', options);

0 10

20 30

40 50

-20 -10 0 10 20 -30 -20 -10 0 10 20 30

提 示 ：

8 若要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下 直接輸入 lorenz 指令即可。非常值得一看！！

一般而言，假設我們的 OutputFcn 設成“myfunc”：

options = odeset('OutputFcn', 'myfunc')

則在積分開始時，ODE 指令會呼叫 myfunc(tspan, y0, 'init')以讓 myfunc 進行各種初始化動作；在每個積分步驟中，ODE 指令將會持續呼叫

三度空 間的相 位圖形

status=myfunc(t, y)，若 status=1，則停止積分；在積分結束時，ODE 指 令會呼叫 myfunc([ ], [ ], 'done')，以讓 myfunc 進行收尾動作。

OutputSel 可指定要傳送到 OutputFcn 的狀態變數之元素。例如，若只 要傳送第一及第三個 Lorenz 渾沌方程式的狀態戀數至 odeplot，可輸 入如下：

>> options = odeset('OutputSel', [1,3]); % 只畫出第一和第三個

% 狀態變數

>> ode45('lorenz1', [0 10], [20 5 -5]', options);

0 2 4 6 8 10

-30 -20 -10 0 10 20 30 40 50

Refine 欄位可使輸出點由內差法而增加數量，以得到較平滑的輸出曲 線。在下例中，我們利用 Refine 欄位來使 ode23 的輸出點個數增為原 先的三倍：

>> subplot(2,1,1); ode23('vdp1', [0 25], [3 3]');

>> options = odeset('Refine', 3);

>> subplot(2,1,2); ode23('vdp1', [0 25], [3 3]', options);

第一個狀態變數

第三個狀態變數

0 5 10 15 20 25 -4

-2 0 2 4

0 5 10 15 20 25

-4 -2 0 2 4

當 Stat=on 時，ODE 指令會在執行完畢後顯示計算過程的各種統計數 字，例如：

>> [t, y] = ode45('vdp1', [0 25], [3 3]', odeset('Stat', 'on'));

71 successful steps 10 failed attempts

487 function evaluations 0 partial derivatives 0 LU decompositions

0 solutions of linear systems

相同的統計數字，也可由 ODE 指令的第三個輸出變數傳回，例如：

>> [t, y, s] = ode45('vdp1', [0 25], [3 3]');

>> s s = 71 10

正常 資料點

資料點 增為 3 偣

487 0 0 0

MaxStep 及 InitialStep 欄位可用來調整最大積分步長（Step Size）及起 始積分步長。但一般而言，您不必去調整這兩個數值，因為 ODE 指 令本身就具有非常先進的步長自動調適功能。特別值得注意的是：

l 若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加。

l 如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟。

其它欄位因較少用到，不再贅述，有興趣的讀者，可直接翻閱 MATLAB 有關 ODE 的手冊。

29-4 ODE 檔案的進階用法

在第二節中，我們已知道基本 ODE 檔案的寫法，如 vdp1.m、vdp2.m 及 lorenz1.m 等。在本節中，將更進一步介紹 ODE 檔案的進階用法，以 使 ODE 指令能夠迅速且準確地算出積分結果。

首先，我們可將 tspan（積分時間範圍）、y0（起始值）及 options（ODE 參數）置於 ODE 檔案中，這些變數必須能由 ODE 檔案傳回，其格式 為：

[tspan, y0, options] = odefile([], [], 'init')

其中我們假設 odefile 即是我們的 ODE 檔案。若 odefile 滿足上述要 求，則我們可以直接呼叫 ODE 指令如下：

[t, y] = solver('odefile')

其中 solver 為前述的任一個 ODE 指令，它可由 odefile 直接得到 tspan、y0 及 options 等積分所需的資訊。

以前述的 van der Pol 為例，若要能夠傳回 tspan、y0 及 options，vdp1.m 須改寫如下（vdp3.m）：

>> type vdp3.m

function [output1, output2, output3] = vdp3(t, y, flag) if strcmp(flag, '')

mu = 1;

output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'),

output1 = [0; 25]; % Time span

output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options

end

此時我們可以直接呼叫 ODE 指令如下：

>> ode45('vdp3');

0 5 10 15 20 25

-4 -2 0 2 4

0 5 10 15 20 25

-4 -2 0 2 4

此外，van der Pol 的微分方程式有一個參數

µ

參數的值，我們可將 vdp3.m 改寫如下（vdp4.m）：

>> type vdp4.m

function [output1, output2, output3] = vdp4(t, y, flag, mu) if nargin < 4 | isempty(mu),

mu = 1;

end

if strcmp(flag, '')

output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'),

output1 = [0; 25]; % Time span

output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options

end

此時

µ

µ

µ

>> subplot(2,1,1); ode45('vdp4', [], [], [], 1); % mu = 1

>> subplot(2,1,2); ode45('vdp4', [], [], [], 3); % mu = 3

0 5 10 15 20 25 -4

-2 0 2 4

0 5 10 15 20 25

-10 -5 0 5 10

在上圖中，

µ

此 外 ， 為 解 決 其 它 較 複 雜 的 ODEs 及 DAEs （ Differential Algebra Equations），ODE 檔案亦可在不同的旗號（Flag）下傳回不同的資訊，

以下是一個完整的列表：

整理：由 odeset 產生的 ODE 選項

旗號 傳回值

（空字串） dy(=F(t,y))

init tspan, yo 及 options

jacobian Jacobian 矩陣 J(t,y) = 2F/2Y jpattern Jacobian sparsity pattern 之矩陣

µ

µ

mass 解 M(t,y)y' = F(t,y) 所須的 質量矩陣 M

events 定義事件發生點的各種資訊

後四者因為不常用到，在此不再贅述。有興趣的讀者可參考 MATLAB 有關 ODE 的手冊。