11-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 結構變數中，某個欄位的值，其格式如下：

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

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

value = odeget(otpions, 'name');

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

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

Example 1: 11-常微分方程式/odeset01.m

odeset AbsTol: [ positive scalar or vector {1e-6} ] RelTol: [ positive scalar {1e-3} ] NormControl: [ on | {off} ] NonNegative: [ vector of integers ] OutputFcn: [ function_handle ] OutputSel: [ vector of integers ] Refine: [ positive integer ] Stats: [ on | {off} ] InitialStep: [ positive scalar ] MaxStep: [ positive scalar ] BDF: [ on | {off} ] MaxOrder: [ 1 | 2 | 3 | 4 | {5} ] Jacobian: [ matrix | function_handle ] JPattern: [ sparse matrix ] Vectorized: [ on | {off} ] Mass: [ matrix | function_handle ] MStateDependence: [ none | {weak} | strong ] MvPattern: [ sparse matrix ] MassSingular: [ yes | no | {maybe} ] InitialSlope: [ vector ] Events: [ function_handle ]

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

類別	欄位名稱	資料型態	預設值	說明
誤差容忍度之相關欄位	RelTol	正純量	$10^{-3}$	相對誤差容忍度
	AbsTo1	正純量或向量	$10^{-6}$	絕對誤差容忍度
積分輸出之相關欄性	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') 傳回 $\partial F/\partial y$，則 Jacobian = 'on'
	Jpattern	on 或 off	off	若 F([ ], [ ], 'JPattern') 傳回 $\partial F/\partial y$，且 $\partial F/\partial y$ 是稀疏矩陣，則 Jpattem = 'on'
	Vectorized	on 或 off	off	若 F(t, [y1, y2…..]) 傳回 [F(t,y1), F(t,y2)…..]，則 Vectorized = 'on'
積分步長（Step Sizes）之相關欄位	Max Step	正純量		ODE 指令之積分步長 的上限
	Initial Step	正純量		起始步長的建議值
質量矩陣之相關欄位	Mass	none，M，M(t)，或 M(t, y)	none	表明 ODE 指令案是否會傳回質量矩陣
	MassSingular	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)，必須滿足下列方程式：

$$ |e(i)| \leq max(RelTol*|y(i)|, AbsTol(i)) $$

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

Example 2: 11-常微分方程式/odeRelTol01.m

subplot(2,1,1); ode45('vdp1', [0 25], [3 3]'); title('RelTol=0.01'); options = odeset('RelTol', 0.00001); subplot(2,1,2); ode45('vdp1', [0 25], [3 3]', options); title('RelTol=0.0001');

在上述範例中，第一個圖所使用的相對誤差值是0.01（預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長。

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

• odephas2：畫出 2-D 的相位平面（Phase Plane）
• odephas3：畫出 3-D 的相位平面
• odeprint：隨時在指令視窗印出計算結果

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

11-常微分方程式/lorenzOde.m

function dy = lorenzOde(t, y) % LORENZODE: 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 對上述Lorenz 渾沌方程式進行數值積分，可見此範例：

Example 3: 11-常微分方程式/odeLorenz01.m

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

請注意在上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形。若要畫三度空間之相位圖形，可使用下列範例：

Example 4: 11-常微分方程式/odeLorenz02.m

options = odeset('OutputFcn', 'odephas3'); % 使用 odephas3 進行繪圖 ode45('lorenzOde', [0 10], [20 5 -5]', options);

上述圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線。

Hint

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

一般而言，假設我們的 OutputFcn 設成"myFunc"：

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

則在整個積分過程中，myFunc 被呼叫的情況如下：

1. 積分開始時，ODE 指令會呼叫 myFunc(tspan, y0, 'init') 以讓 myFunc 進行各種初始化動作。
2. 在每個積分步驟中，ODE 指令將會持續呼叫 status=myFunc(t, y)，若 status=1，則停止積分。
3. 在積分結束時，ODE 指令會呼叫 myFunc([], [], 'done')，以讓 myFunc 進行收尾動作。

因此，若我們要在積分過程中產生動畫，就可以根據上述原則來設計 myFunc.m。

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

Example 5: 11-常微分方程式/odeOutputSelect01.m

options = odeset('OutputSel', [1,3]); % 只畫出第一和第三個狀態變數 ode45('lorenzOde', [0 10], [20 5 -5]', options);

在上圖中，我們只畫出了第一和第三個狀態變數。

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

Example 6: 11-常微分方程式/odeRefine01.m

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

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

Example 7: 11-常微分方程式/odeShowStats01.m

[t, y] = ode45('vdp1', [0 25], [3 3]', odeset('Stat', 'on'));71 successful steps 10 failed attempts 487 function evaluations

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

Example 8: 11-常微分方程式/odeShowStats02.m

[t, y, s] = ode45('vdp1', [0 25], [3 3]'); ss = 71 10 487 0 0 0

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

• 若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加。
• 如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟。

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

---

MATLAB程式設計：進階篇