MATLAB 程式設計進階篇常微分方程式.ppt

《MATLAB 程式設計進階篇常微分方程式.ppt》由会员分享，可在线阅读，更多相关《MATLAB 程式設計進階篇常微分方程式.ppt（46页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、MATLAB 程式設計進階篇常微分方程式,張智星 jangcs.nthu.edu.tw http:/www.cs.nthu.edu.tw/jang 清大資工系 多媒體檢索實驗室,11-1 ODE 指令列表,MATLAB 用於求解常微分方程式的指令:,ODE 指令列表,指令項目繁多，最主要可分兩大類 適用於 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接選用 ode45、ode23 或 ode113 來求解 適用Stiff 系統 速率（即微分值）差異相常大 使用一般的 ode45、ode23 或 ode113 來求解，可能會使得積分的步長（Step Sizes）變得

2、很小，以便降低積分誤差至可容忍範圍以內，會導致計算時間過長 專門對付 Stiff 系統的指令，例如 ode15s、ode23s、ode23t 及 ode 23tb,提示,使用 Simulink 來求解常微分方程式 Simulink是和MATLAB共同使用的一套軟體 可使用拖拉的方式來建立動態系統 可直接產生C程式碼或進行動畫顯示 功能非常強大,11-2 ODE 指令基本用法,使用 ODE 指令時，必須先將要求解的 ODE 表示成一個函式 輸入為 t（時間）及 y（狀態變數，State Variables） 輸出則為 dy（狀態變數的微分值） ODE 函式的檔名為 odeFile.m，則呼叫 O

3、DE 指令的格式如下： t, y = solver (odeFile, t0, t1, y0); t0, t1 是積分的時間區間 y0 代表起始條件（Initial Conditions） solver 是前表所列的各種 ODE 指令 t 是輸出的時間向量 y 則是對應的狀態變數向量,ODE 指令基本用法,以 van der Pol 微分方程為例，其方程式為： 化成標準格式 可用向量來表示成一般化的數學式 為一向量，代表狀態變數,ODE 指令基本用法,假設 =1， ODE 檔案（vdp1.m）可顯示如下： type vdp1.m function dy = vdp1(t, y) mu = 1;

4、 dy = y(2); mu*(1-y(1)2)*y(2)-y(1); 有了 ODE 檔案，即可選用前述之 ODE 指令來求解,ODE 指令基本用法範例-1 (I),在 =1 時，van der Pol 方程式並非 Stiff 系統，所以使用ode45來畫出積分的結果 範例11-1：odeBasic01.m,ode45(vdp1, 0 25, 3 3);,ODE 指令基本用法範例-1 (II),0, 25 代表積分的時間區間，3 3 則代表起始條件（必須以行向量來表示） 因為沒有輸出變數，所以上述程式執行結束後，MATLAB 只會畫出狀態變數對時間的圖形,ODE 指令基本用法範例-2 (I),

5、要取得積分所得的狀態變數及對應的時間，可以加上輸出變數以取得這些資料 範例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);,ODE 指令基本用法範例-2 (II),ODE 指令基本用法範例-3 (I),也可以畫出 及 在 相位平面（Phase Plane ）的運動情況 範例11-3：odePhastPlot01.m,t, y = ode45(vdp1,

6、 0 25, 3 3); plot(y(:,1), y(:,2), -o); xlabel(y(t); ylabel(y(t);,ODE 指令基本用法範例-3 (II),當 值越來越大時，van der Pol 方程式就漸漸變成一個 Stiff 的系統，此時若要解此方程式，就必須改用專門對付 Stiff 系統的指令,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);,ODE 指令基本

7、用法範例-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);% 注意單引號的重覆使用,ODE 指令基本用法範例-4 (III),由上圖可知， 的變化相當劇烈（超過 ），這就是 St

8、iff 系統的特色,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),ODE 指令基本用法範例-5 (II),

9、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，依此類推 未被設定的欄位，其欄位值即為預設值,ODE

10、 指令的選項,也可以只改變一個現存的 options 結構變數中，某個欄位的值，其格式如下： newOptions = odeset(options, name, value); 若要讀出某個欄位的值，可用 odeget，其格式如下： value = odeget(otpions, name); 其中 name 為欄位名稱，value 即為對應的欄位值 當使用 odeset 指令時，若無任何輸入變數，則 odeset 會顯示所有的欄位名稱及欄位值，並以大括號代表預設值,ODE 指令的選項, odeset AbsTol: positive scalar or vector 1e-6 RelTol

11、: 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_h

12、andle 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 ,由 odeset 產生的 ODE 選項,由 odeset 產生的 ODE 選項,若F(t, y, Jacobian) 傳回,若 F( , , JPatte

13、rn) 傳回,，且,由 odeset 產生的 ODE 選項,常用到的欄位來進行說明,f 在積分誤差容忍度方面，每一次積分所產生的局部誤差 e(i)，必須滿足下列方程式： max(RelTol* , AbsTol(i) 其中i 代表第 i 個狀態變數 降低 RelTol 及 AbsTol 來求得更精確的積分結果,範例-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); ode

14、45(vdp1, 0 25, 3 3, options); title(RelTol=0.0001);,範例-1 (II),第一個圖所使用的相對誤差值是0.01（預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長,積分輸出方面說明,積分輸出的相關處理方面 選用一個 OutputFcn 當 ODE 指令沒有輸出變數時，此輸出函式 OutputFcn 會被 MATLAB 呼叫 OutputFcn 的預設值是”odeplot”，其功能為畫出所有的狀態變數 其它可用的函式 odephas2：畫出 2-D 的相位平面（Phase Plane） ode

15、phas3：畫出 3-D 的相位平面 odeprint：隨時在指令視窗印出計算結果,積分輸出方面說明,以 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;,範例-2,使用 ode45 對Lo

16、renz 渾沌方程式進行數值積分 範例11-7：odeLorenz01.m 上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形,ode45(lorenzOde, 0 10, 20 5 -5);,範例-3,上述範例畫三度空間之相位圖形 範例11-8：odeLorenz02.m 圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線,options = odeset(OutputFcn, odephas3); % 使用 odephas3 進行繪圖 ode45(lorenzOde, 0 10, 20 5 -5, options);,提示,要觀看 Lorenz

17、渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入lorenz 指令,範例-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 的狀

18、態變數之元素,範例-4 (II),只要傳送第一及第三個 Lorenz 渾沌方程式的狀態變數至 odeplot 範例11-9：odeOutputSelect01.m,options = odeset(OutputSel, 1,3) % 只畫出第一和第三個狀態變數 ode45(lorenzOde, 0 10, 20 5 -5, options);,範例-5 (I),Refine 欄位可以使用內差法來增加輸出狀態變數的密度，以得到較平滑的輸出曲線 用 Refine 欄位使 ode23 的輸出點個數增為原先三倍： 範例11-10：odeRefine01.m,subplot(2,1,1); ode23(

19、vdp1, 0 25, 3 3); subplot(2,1,2); options = odeset(Refine, 3); ode23(vdp1, 0 25, 3 3, options);,範例-5 (II),範例-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);,範例-7,相同的統計數字，

20、也可由 ODE 指令的第三個輸出變數傳回 範例11-12：odeShowStats02.m s = 71 10 487 0 0 0,t, y, s = ode45(vdp1, 0 25, 3 3); s,說明,MaxStep 及 InitialStep 欄位可用來調整最大積分步長及起始積分步長 一般而言，不必去調整這兩個數值，因為 ODE 指令本身就具有步長自動調適功能 若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加 如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降

21、MaxStep 是最後的步驟,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 指令，它可由 od

22、eFile 直接得到 tspan、y0 及 options 等積分所需的資訊,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

23、 = 0; 25;% Time span output2 = 3; 3;% Initial conditions output3 = odeset;% ODE options end,ODE 檔案的進階用法範例-1 (II),範例11-13：odeAdvanced01.m,ode45(vdp3),ODE 檔案的進階用法範例-2 (I),van der Pol 的微分方程式有一個參數 ，希望從外面傳入此參數的值(vdp4.m ) type vdp4.m function output1, output2, output3 = vdp4(t, y, flag, mu) if nargin 4 |

24、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,ODE 檔案的進階用法範例-2 (II),就變成一個選擇性（Optional）的參數，其預設值為 1 將 的值從 MATLAB 傳入，並畫出不同 值下的 van der

25、 Pol 方程式的狀態變數： 範例11-14：odeAdvanced02.m,subplot(2,1,1); ode45(vdp4, , , , 1);% mu=1 subplot(2,1,2); ode45(vdp4, , , , 3);% mu=3,ODE 檔案的進階用法範例-2 (III),的值分別是 1 及 3 用到了許多空矩陣，這些空矩陣代表取用預設值，因此 ode45 會直接從 vdp4.m 取用時間區間及變數起始值 也可以傳入二個或更多的參數，MATLAB 及 ODE 指令對於可傳入的參數個數並無設限,ODE 檔案的進階用法列表,為解決其它較複雜的 ODEs 及 DAEs（Differential Algebra Equations），ODE 檔案亦可在不同的旗號（Flag）下傳回不同的資訊，列表如下：,