MATLAB 程式設計進階篇 常微分方程式 張智星 清大資工系 多媒體檢索實驗室.

Presentation on theme: "MATLAB 程式設計進階篇 常微分方程式 張智星 清大資工系 多媒體檢索實驗室."— Presentation transcript:

1 MATLAB 程式設計進階篇 常微分方程式 張智星 jang@cs.nthu.edu.tw http://www.cs.nthu.edu.tw/~jang 清大資工系 多媒體檢索實驗室

2 MATLAB 程式設計進階篇：常微分方程 式 11-1 ODE 指令列表 MATLAB 用於求解常微分方程式的指令 : 指令方法適用 ODE 類別 ode45Explicit Runge-Kutta (4, 5) pair of Dormand-PrinceNonstiff ODE ode23Explicit Runge-Kutta (2, 3) pair of Bogacki and ShampineNonstiff ODE ode113Variable order Adams-Bashforth-Moulton PECE solverNonstiff ODE ode15sNumerical differentiation formulas (NDFS)Stiff ODE ode23sModified Rosenbrock formula of order 2Stiff ODE ode23tTrapezoidal rule with a “free” interpolantStiff ODE ode23tb Implicit Runge-Kutta formula with a backward differentiation formula of order two Stiff ODE

3 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令列表 指令項目繁多，最主要可分兩大類 適用於 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接選用 ode45 、 ode23 或 ode113 來求解 適用 Stiff 系統 速率（即微分值）差異相常大 使用一般的 ode45 、 ode23 或 ode113 來求解，可能會使 得積分的步長（ Step Sizes ）變得很小，以便降低積分誤差 至可容忍範圍以內，會導致計算時間過長 專門對付 Stiff 系統的指令，例如 ode15s 、 ode23s 、 ode23t 及 ode 23tb

4 MATLAB 程式設計進階篇：常微分方程 式 提示 使用 Simulink 來求解常微分方程式 Simulink 是和 MATLAB 共同使用的一套軟體 可使用拖拉的方式來建立動態系統 可直接產生 C 程式碼或進行動畫顯示 功能非常強大

5 MATLAB 程式設計進階篇：常微分方程 式 11-2 ODE 指令基本用法 使用 ODE 指令時，必須先將要求解的 ODE 表示成 一個函式 輸入為 t （時間）及 y （狀態變數， State Variables ） 輸出則為 dy （狀態變數的微分值） ODE 函式的檔名為 odeFile.m ，則呼叫 ODE 指令的 格式如下： [t, y] = solver ('odeFile', [t0, t1], y0); [t0, t1] 是積分的時間區間 y0 代表起始條件（ Initial Conditions ） solver 是前表所列的各種 ODE 指令 t 是輸出的時間向量 y 則是對應的狀態變數向量

6 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法 以 van der Pol 微分方程為例，其方程式為： 化成標準格式 可用向量來表示成一般化的數學式 為一向量，代表狀態變數

7 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法 假設 =1 ， ODE 檔案（ vdp1.m ）可顯示如下： >> type vdp1.m function dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; 有了 ODE 檔案，即可選用前述之 ODE 指令來求解

8 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -1 (I) 在 =1 時， van der Pol 方程式並非 Stiff 系統， 所以使用 ode45 來畫出積分的結果 範例 11-1 ： odeBasic01.m ode45('vdp1', [0 25], [3 3]');

9 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -1 (II) [0, 25] 代表積分的時 間區間， [3 3] ’ 則代表 起始條件（必須以行向 量來表示） 因為沒有輸出變數，所 以上述程式執行結束後， MATLAB 只會畫出狀態 變數對時間的圖形

10 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -2 (I) 要取得積分所得的狀態變數及對應的時間，可以加 上輸出變數以取得這些資料 範例 11-2 ： 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)');

11 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -2 (II)

12 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -3 (I) 也可以畫出 及 在 相位平面（ Phase Plane ） 的運動情況 範例 11-3 ： odePhastPlot01.m [t, y] = ode45('vdp1', [0 25], [3 3]'); plot(y(:,1), y(:,2), '-o'); xlabel('y(t)'); ylabel('y''(t)');

13 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -3 (II) 當 值越來越大時， van der Pol 方程式就漸 漸變成一個 Stiff 的系 統，此時若要解此方程 式，就必須改用專門對 付 Stiff 系統的指令

14 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -4 (I) 將 值改成 1000 ， ODE 檔案改寫成（ vdp2.m ）： >> type vdp2.m function dy = vdp2(t, y) mu = 1000; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

15 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -4 (II) 選用專門對付 Stiff 系統的 ODE 指令，例如 ode15s ，來求解此系統並作圖顯示 範例 11-4 ： 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)'); % 注意單引號「 ' 」的重覆使用

16 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -4 (III) 由上圖可知， 的變化相當劇烈（超過 ）， 這就是 Stiff 系統的特色

17 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -5 (I) 若要畫出二維平面相位圖，可以使用下列範例： 範例 11-5 ： ode15s02.m 若要產生在某些特定時間點的狀態變數值，則呼叫 ODE 指令的格式可改成： [t, y] = solver('odeFile', [t0, t1, …, tn], y0); 其中 [t0, t1, …, tn] 即是特定時間點所形成的向量 [t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(1,1,1); plot(y(:, 1), y(:, 2), '-o'); xlabel('y(t)'); ylabel('y''(t)')

18 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令基本用法範例 -5 (II)

19 MATLAB 程式設計進階篇：常微分方程 式 11-3 ODE 指令的選項 ODE 指令可以接受第四個輸入變數，代表積分過程 用到的各種選項（ Options ），此種 ODE 指令的格 式為： [t, y] = solver('odeFile', [t0, tn], y0, options); 其中 options 是由 odeset 指令來控制的結構變數 結構變數即包含了積分過程用到的各種選項 odeset 的一般格式如下： options = odeset('name1', value1, 'name2', value2, … ) 其欄位 name1 的值為 value1 ，欄位 name2 的值為 value2 ， 依此類推 未被設定的欄位，其欄位值即為預設值

20 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令的選項 也可以只改變一個現存的 options 結構變數中，某 個欄位的值，其格式如下： newOptions = odeset(options, 'name', value); 若要讀出某個欄位的值，可用 odeget ，其格式如下： value = odeget(otpions, 'name'); 其中 name 為欄位名稱， value 即為對應的欄位值 當使用 odeset 指令時，若無任何輸入變數，則 odeset 會顯示所有的欄位名稱及欄位值，並以大括 號代表預設值

21 MATLAB 程式設計進階篇：常微分方程 式 ODE 指令的選項 >> 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 ]

22 MATLAB 程式設計進階篇：常微分方程 式 由 odeset 產生的 ODE 選項 類別欄位名稱資料型態預設值說明 誤差容忍度之相 關欄位 RelTol 正純量相對誤差容忍度 AbsTo1 正純量或向量絕對誤差容忍度 積分輸出之相關 欄性 OutPutFcn 字串 ‘odeplot’ 輸出函式（若 ODE 指令無 輸出變數，則在數值積分執 行完畢後， MATLAB 會呼叫 此輸出函式） OutputSel 索引向量全部 ODE 指令之輸出變數的索引 值，以決定那些輸出變數之 元素將被送到輸出函式 Refine 正整數 1 或 4 (for ode45) Refine = 2 可產生兩倍數量 的輸出點， Refine = 3 可產 生三倍數量的輸出點，依此 類推。 Stats on 或 off off Stats = 'on' 產生計算過程的 各種統計資料。

23 MATLAB 程式設計進階篇：常微分方程 式 由 odeset 產生的 ODE 選項 若 F(t, y, ’ Jacobian') 傳回 若 F([ ], [ ], ’ JPattern ’ ) 傳回 ，且，且 Jacobian 矩陣 之相關欄位 Jconstant on 或 off off 如果 Jacobian 矩陣常數， 則 JConstant = 'on' Jacobian on 或 off off ，則 Jacobian = 'on‘ Jpattern on 或 off off 是稀疏矩陣， 則 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 正純量起始步長的建議值

24 MATLAB 程式設計進階篇：常微分方程 式 由 odeset 產生的 ODE 選項 質量矩陣之相關 欄位 Mass none ， M ， M(t) ，或 M(t, y) none 表明 ODE 指令案是否會傳 回質量矩陣 MassSing ular yes ， no 或 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'

25 MATLAB 程式設計進階篇：常微分方程 式 常用到的欄位來進行說明 f 在積分誤差容忍度方面，每一次積分所產生的局 部誤差 e(i) ，必須滿足下列方程式： max(RelTol*, AbsTol(i)) 其中 i 代表第 i 個狀態變數 降低 RelTol 及 AbsTol 來求得更精確的積分結果

26 MATLAB 程式設計進階篇：常微分方程 式 範例 -1 (I) 範例 11-6 ： 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');

27 MATLAB 程式設計進階篇：常微分方程 式 範例 -1 (II) 第一個圖所使用的相對 誤差值是 0.01 （預設 值），第二個圖所使用 的相對誤差值是 0.00001 ，因此我們得 到較細密的點，但所花 的計算時間也會比較長

28 MATLAB 程式設計進階篇：常微分方程 式 積分輸出方面說明 積分輸出的相關處理方面 選用一個 OutputFcn 當 ODE 指令沒有輸出變數時，此輸出函式 OutputFcn 會 被 MATLAB 呼叫 OutputFcn 的預設值是 ” odeplot ” ，其功能為畫出所有的狀 態變數 其它可用的函式 odephas2 ：畫出 2-D 的相位平面（ Phase Plane ） odephas3 ：畫出 3-D 的相位平面 odeprint ：隨時在指令視窗印出計算結果

29 MATLAB 程式設計進階篇：常微分方程 式 積分輸出方面說明 以 Lorenz 渾沌方程式（ Lorenz Chaotic Equation ） 為例 >> type lorenzOde.m function dy = lorenzOde(t, y) % LORENZODE: ODE file for Lorenz chaotic equation IGMA = 10.; RHO = 28; BETA = 8./3.; A = [ -BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1 ]; dy = A*y;

30 MATLAB 程式設計進階篇：常微分方程 式 範例 -2 使用 ode45 對 Lorenz 渾沌方程式進行數值積分 範例 11-7 ： odeLorenz01.m 上列圖中，共有三條曲線，代表三個狀態變數隨 時間變化的圖形 ode45('lorenzOde', [0 10], [20 5 -5]');

31 MATLAB 程式設計進階篇：常微分方程 式 範例 -3 上述範例畫三度空間之相位圖形 範例 11-8 ： odeLorenz02.m 圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、 以時間為參數的一條三度空間中的曲線 options = odeset('OutputFcn', 'odephas3'); % 使用 odephas3 進行繪圖 ode45('lorenzOde', [0 10], [20 5 -5]', options);

32 MATLAB 程式設計進階篇：常微分方程 式 提示 要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入 lorenz 指令

33 MATLAB 程式設計進階篇：常微分方程 式 範例 -4 (I) 假設 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 的狀態變數之元素

34 MATLAB 程式設計進階篇：常微分方程 式 範例 -4 (II) 只要傳送第一及第三個 Lorenz 渾沌方程式的狀態變 數至 odeplot 範例 11-9 ： odeOutputSelect01.m options = odeset('OutputSel', [1,3]) % 只畫出第一和第三個狀態變數 ode45('lorenzOde', [0 10], [20 5 -5]', options);

35 MATLAB 程式設計進階篇：常微分方程 式 範例 -5 (I) Refine 欄位可以使用內差法來增加輸出狀態變數的 密度，以得到較平滑的輸出曲線 用 Refine 欄位使 ode23 的輸出點個數增為原先三倍： 範例 11-10 ： 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);

36 MATLAB 程式設計進階篇：常微分方程 式 範例 -5 (II)

37 MATLAB 程式設計進階篇：常微分方程 式 範例 -6 當 Stat=on 時， ODE 指令會在執行完畢後顯示計 算過程的各種統計數字 範例 11-11 ： odeShowStats01.m 71 successful steps 10 failed attempts 487 function evaluations [t, y] = ode45('vdp1', [0 25], [3 3]', odeset('Stat', 'on'));

38 MATLAB 程式設計進階篇：常微分方程 式 範例 -7 相同的統計數字，也可由 ODE 指令的第三個輸出 變數傳回 範例 11-12 ： odeShowStats02.m s = 71 10 487 0 [t, y, s] = ode45('vdp1', [0 25], [3 3]'); s

39 MATLAB 程式設計進階篇：常微分方程 式 說明 MaxStep 及 InitialStep 欄位可用來調整最大積分步 長及起始積分步長 一般而言，不必去調整這兩個數值，因為 ODE 指令本身 就具有步長自動調適功能 若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大 幅增加 如果積分結果不甚準確，請勿先調降 MaxStep ，您應先調 降 RelTol 及 AbsTol 。調降 MaxStep 是最後的步驟

40 MATLAB 程式設計進階篇：常微分方程 式 11-4 ODE 檔案的進階用法 更進一步介紹 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 等積分所需的資訊

41 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法範例 -1 (I) 以前述的 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

42 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法範例 -1 (II) 範例 11-13 ： odeAdvanced01.m ode45('vdp3')

43 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法範例 -2 (I) van der Pol 的微分方程式有一個參數 ，希望從外 面傳入此參數的值 (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

44 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法範例 -2 (II) 就變成一個選擇性（ Optional ）的參數，其預設 值為 1 將 的值從 MATLAB 傳入，並畫出不同 值下的 van der Pol 方程式的狀態變數： 範例 11-14 ： odeAdvanced02.m subplot(2,1,1); ode45('vdp4', [], [], [], 1);% mu=1 subplot(2,1,2); ode45('vdp4', [], [], [], 3);% mu=3

45 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法範例 -2 (III) 的值分別是 1 及 3 用到了許多空矩陣，這 些空矩陣代表「取用預 設值」，因此 ode45 會 直接從 vdp4.m 取用時 間區間及變數起始值 也可以傳入二個或更多 的參數， MATLAB 及 ODE 指令對於可傳入的 參數個數並無設限

46 MATLAB 程式設計進階篇：常微分方程 式 ODE 檔案的進階用法列表 為解決其它較複雜的 ODEs 及 DAEs （ Differential Algebra Equations ）， ODE 檔案亦可在不同的旗號 （ Flag ）下傳回不同的資訊，列表如下： 旗號傳回值 （空字串） 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 定義事件發生點的各種資訊