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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 程式設計進階篇常微分方程式 程式 設計進階篇常 微分 方程式
限制150内