所有的 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 會顯示所有的欄位名稱及欄位值,並以大括號代表預設值,例如:
這些欄位名稱及相關說明,可見下表:
類別 欄位名稱 資料型態 預設值 說明 誤差容忍度之相關欄位 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 來求得更精確的積分結果。例如:
在上述範例中,第一個圖所使用的相對誤差值是0.01(預設值),第二個圖所使用的相對誤差值是0.00001,因此我們得到較細密的點,但所花的計算時間也會比較長。
在積分輸出的相關處理方面,我們可以選用一個 OutputFcn,使得當 ODE 指令沒有輸出變數時,此輸出函式 OutputFcn 會被 MATLAB 呼叫。OutputFcn 的預設值是”odeplot”,其功能為畫出所有的狀態變數,其它可用的函式還有:
- odephas2:畫出 2-D 的相位平面(Phase Plane)
- odephas3:畫出 3-D 的相位平面
- odeprint:隨時在指令視窗印出計算結果
以 Lorenz 渾沌方程式(Lorenz Chaotic Equation)為例,其 ODE 檔案之內容可顯示如下:
欲使用 ode45 對上述Lorenz 渾沌方程式進行數值積分,可見此範例:
請注意在上列圖中,共有三條曲線,代表三個狀態變數隨時間變化的圖形。若要畫三度空間之相位圖形,可使用下列範例:
上述圖形中只出現一條曲線,此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線。
一般而言,假設我們的 OutputFcn 設成"myFunc":
options = odeset('OutputFcn', 'myFunc') 則在整個積分過程中,myFunc 被呼叫的情況如下:
因此,若我們要在積分過程中產生動畫,就可以根據上述原則來設計 myFunc.m。
- 積分開始時,ODE 指令會呼叫 myFunc(tspan, y0, 'init') 以讓 myFunc 進行各種初始化動作。
- 在每個積分步驟中,ODE 指令將會持續呼叫 status=myFunc(t, y),若 status=1,則停止積分。
- 在積分結束時,ODE 指令會呼叫 myFunc([], [], 'done'),以讓 myFunc 進行收尾動作。
OutputSel 可指定要傳送到 OutputFcn 的狀態變數之元素。例如,若只要傳送第一及第三個 Lorenz 渾沌方程式的狀態戀數至 odeplot,可輸入如下:
在上圖中,我們只畫出了第一和第三個狀態變數。
Refine 欄位可以使用內插法來增加輸出狀態變數的密度,以得到較平滑的輸出曲線。在下例中,我們利用 Refine 欄位來使 ode23 的輸出點個數增為原先的三倍:
當 Stat=on 時,ODE 指令會在執行完畢後顯示計算過程的各種統計數字,例如:
相同的統計數字,也可由 ODE 指令的第三個輸出變數傳回,例如:
MaxStep 及 InitialStep 欄位可用來調整最大積分步長及起始積分步長。但一般而言,您不必去調整這兩個數值,因為 ODE 指令本身就具有步長自動調適功能。特別值得注意的是:
- 若要產生更多輸出點,可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果,但是計算時間可能會大幅增加。
- 如果積分結果不甚準確,請勿先調降 MaxStep,您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟。
其它欄位因較少用到,不再贅述,有興趣的讀者,可直接翻閱 MATLAB 有關 ODE 的手冊。
MATLAB程式設計:進階篇