常微分方程式ppt课件

1、MATLAB 程序设计进阶篇常微分方程式張智星.tw/jang清大资工系 多媒体检索实验室MATLAB 程式設計進階篇：常微分方程式11-1 ODE 指令列表 nMATLAB用于求解常微分方程式的指令: 指令方法適用ODE類別ode45Explicit Runge-Kutta (4, 5) pair of Dormand-PrinceNonstiff ODEode23Explicit Runge-Kutta (2, 3) pair of Bogacki and ShampineNonstiff ODEode113Variable ord

2、er Adams-Bashforth-Moulton PECE solverNonstiff ODEode15sNumerical differentiation formulas (NDFS)Stiff ODEode23sModified Rosenbrock formula of order 2Stiff ODEode23t Trapezoidal rule with a “free” interpolantStiff ODEode23tbImplicit Runge-Kutta formula with a backward differentiation formula of orde

3、r twoStiff ODEMATLAB 程式設計進階篇：常微分方程式ODE 指令列表n指令项目繁多，最主要可分两大类n适用于Nonstiff系统n一般的常微分方程式都是Nonstiff系统n直接选用ode45、ode23或ode113来求解n适用Stiff系统n速率即微分值差异相常大n使用一般的ode45、ode23或ode113来求解，可能会使得积分的步长Step Sizes变得很小，以便降低积分误差至可容忍范围以内，会导致计算时间过长n专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode 23tbMATLAB 程式設計進階篇：常微分方程式提示 n使用 Sim

4、ulink 來求解常微分方程式nSimulink是和MATLAB共同使用的一套软件n可使用拖拉的方式来建立动态系统n可直接产生C程序代码或进行动画显示n功能非常强大 MATLAB 程式設計進階篇：常微分方程式11-2 ODE 指令基本用法 n使用ODE指令时，必须先将要求解的ODE表示成一个函式n输入为t时间及y状态变量，State Variables）n输出则为dy状态变量的微分值）nODE函式的档名为odeFile.m，则呼叫ODE指令的格式如下：nt, y = solver (odeFile, t0, t1, y0);nt0, t1是积分的时间区间ny0代表起始条件Initial Con

5、ditions）nsolver是前表所列的各种ODE指令nt是输出的时间向量ny则是对应的状态变量向量MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法n以van der Pol微分方程为例，其方程式为：n化成标准格式n可用向量来表示成一般化的数学式n 为一向量，代表状态变量0)1 ( 2yyyy1221221)1 (yyyyyy),(tyFyyMATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法n假设=1，ODE档案vdp1.m可显示如下：ntype vdp1.m n function dy = vdp1(t, y)n mu = 1;n dy = y(2); mu*(1

6、-y(1)2)*y(2)-y(1);n有了ODE档案，即可选用前述之ODE指令来求解MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-1 (I)n在 =1 時，van der Pol 方程式并非Stiff系统，所以使用ode45来画出积分的结果n範例11-1：odeBasic01.mode45(vdp1, 0 25, 3 3);MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-1 (II)n 0, 25 代表积分的时间区间，3 3 则代表起始条件必须以行向量来表示）n因为没有输出变量，所以上述程序执行结束后，MATLAB只会画出状态变量对时间的图形05101

7、52025-3-2-101234MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-2 (I)n要取得积分所得的状态变量及对应的时间，可以加上输出变量以取得这些数据n范例11-2：odeGetData01.mt, 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);MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-2 (II)0510152025-3-2-101234Tim

8、e tSolution y(t) and y(t) y(t)y(t)MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-3 (I)n也可以画出 及 在 相位平面Phase Plane的运动情况n範例11-3：odePhastPlot01.m)(ty)( tyt, y = ode45(vdp1, 0 25, 3 3);plot(y(:,1), y(:,2), -o);xlabel(y(t); ylabel(y(t);MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-3 (II)n当 值越来越大时，van der Pol方程式就渐渐变成一个Stiff的系统，此时若

9、要解此方程式，就必须改用专门对付Stiff系统的指令-3-2-101234-3-2-10123y(t)y(t)MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-4 (I)n将 值改成1000，ODE档案改写成vdp2.m）：n type vdp2.m n function dy = vdp2(t, y)n mu = 1000;n dy = y(2); mu*(1-y(1)2)*y(2)-y(1); MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-4 (II)n选用专门对付Stiff系统的ODE指令，例如ode15s，来求解此系统并作图显示n範例11-4：o

10、de15s01.mt, 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);% 注意單引號的重覆使用MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-4 (III)n由上图可知， 的变化相当剧烈超越 ），这就是Stiff系统的特色050010001500200025003000-4-2024Time ty(t)

11、050010001500200025003000-2000-1000010002000Time ty(t)( ty1000MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-5 (I)n若要画出二维平面相位图，可以使用下列范例：n范例11-5：ode15s02.mn若要产生在某些特定时间点的状态变量值，则呼叫ODE指令的格式可改成：n t, y = solver(odeFile, t0, t1, tn, y0);n其中t0, t1, tn即是特定时间点所形成的向量t, y= ode15s(vdp2, 0 3000, 2 1);subplot(1,1,1);plot(y(:, 1

12、), y(:, 2), -o);xlabel(y(t); ylabel(y(t)MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-5 (II)-2.5-2-1.5-1-0.500.511.522.5-1500-1000-500050010001500y(t)y(t)MATLAB 程式設計進階篇：常微分方程式11-3 ODE 指令的選項nODE指令可以接受第四个输入变量，代表积分过程用到的各种选项Options），此种ODE指令的格式为：nt, y = solver(odeFile, t0, tn, y0, options);n其中options是由odeset指令来控制的结构

13、变量n结构变量即包含了积分过程用到的各种选项nodeset的一般格式如下：n options = odeset(name1, value1, name2, value2,)n其字段name1的值为value1，字段name2的值为value2，依此类推n未被设定的字段，其域值即为默认值MATLAB 程式設計進階篇：常微分方程式ODE 指令的選項n也可以只改变一个现存的options结构变量中，某个字段的值，其格式如下：n newOptions = odeset(options, name, value);n若要读出某个字段的值，可用odeget，其格式如下：nvalue = odeget(ot

14、pions, name);n其中name为域名，value即为对应的域值n当使用odeset指令时，若无任何输入变量，则odeset会显示所有的域名及域值，并以大括号代表默认值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 inte

15、gers 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 MvPatter

16、n: sparse matrix MassSingular: yes | no | maybe InitialSlope: vector Events: function_handle MATLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項 310610類別欄位名稱資料型態預設值說明誤差容忍度之相關欄位RelTol正純量相對誤差容忍度AbsTo1正純量或向量絕對誤差容忍度積分輸出之相關欄性OutPutFcn字串odeplot輸出函式（若 ODE 指令無輸出變數，則在數值積分執行完畢後，MATLAB 會呼叫此輸出函式）OutputSel索引向量全部ODE 指令之輸出變數

17、的索引值，以決定那些輸出變數之元素將被送到輸出函式Refine正整數1或4 (for ode45)Refine = 2 可產生兩倍數量的輸出點，Refine = 3 可產生三倍數量的輸出點，依此類推。Statson 或 offoffStats = on 產生計算過程的各種統計資料。MATLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項若F(t, y, Jacobian) 傳回 yF假设 F( , , JPattern) 傳回yF，且yFJacobian 矩陣之相關欄位Jconstanton 或 offoff如果 Jacobian 矩陣常數，則 JConstant =

18、onJacobianon 或 offoff，則 Jacobian = onJpatternon 或 offoff 是稀疏矩陣，則 JPattem = onVectorizedon 或 offoff若 F(t, y1, y2.) 傳回 F(t,y1), F(t,y2).，則 Vectorized = on積分步長（Step Sizes）之相關欄位Max Step正純量ODE 指令之積分步長 的上限Initial Step正純量起始步長的建議值MATLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項質量矩陣之相關欄位Massnone，M，M(t)，或 M(t, y)none

19、表明 ODE 指令案是否會傳回質量矩陣MassSingular yes，no 或 maybemaybe表明質量矩陣是否為 Singular事件發生時間之相關欄位Eventson 或 offoff若 ODE 檔案並傳回事件函式及相關資訊，則 Events = onOde15s 之相關欄位MaxOrder1, 2, 3, 4 或 55積分公式用到的最高次數BDFon 或 offoff若使用 BDF（Backward Differentiation Formula）則 BDF = onMATLAB 程式設計進階篇：常微分方程式常用到的欄位來進行說明 nf在积分误差容忍度方面，每一次积分所产生的局部误

20、差e(i)，必须满足下列方程式：nmax(RelTol* , AbsTol(i)n其中i代表第i个状态变量n降低RelTol及AbsTol来求得更精确的积分结果)(ie)(iyMATLAB 程式設計進階篇：常微分方程式範例-1 (I)n範例11-6：odeRelTol01.msubplot(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);MAT

21、LAB 程式設計進階篇：常微分方程式範例-1 (II)n第一個圖所使用的相對誤差值是0.01預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長 0510152025-4-2024RelTol=0.010510152025-4-2024RelTol=0.0001MATLAB 程式設計進階篇：常微分方程式積分輸出方面說明n积分输出的相关处理方面n选用一个OutputFcnn当ODE指令没有输出变量时，此输出函式OutputFcn会被MATLAB呼叫nOutputFcn的默认值是odeplot”，其功能为画出所有的状态变量n其它可用的函式node

22、phas2：画出2-D的相位平面Phase Plane）nodephas3：画出3-D的相位平面nodeprint：随时在指令窗口印出计算结果MATLAB 程式設計進階篇：常微分方程式積分輸出方面說明n以 Lorenz 渾沌方程式Lorenz Chaotic Equation為例 n type lorenzOde.m n function dy = lorenzOde(t, y)n % LORENZODE: ODE file for Lorenz chaotic equationn IGMA = 10.;n RHO = 28;n BETA = 8./3.;n A = -BETA 0 y(2)n

23、 0 -SIGMA SIGMAn -y(2) RHO -1 ;n dy = A*y;MATLAB 程式設計進階篇：常微分方程式範例-2 n使用 ode45 對Lorenz 渾沌方程式進行數值積分n範例11-7：odeLorenz01.mn上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形 ode45(lorenzOde, 0 10, 20 5 -5);012345678910-30-20-1001020304050MATLAB 程式設計進階篇：常微分方程式範例-3n上述範例畫三度空間之相位圖形 n範例11-8：odeLorenz02.m n圖形中只出現一條曲線，此曲線代表以三個狀態變數為

24、座標、以時間為參數的一條三度空間中的曲線 options = odeset(OutputFcn, odephas3); % 使用 odephas3 進行繪圖ode45(lorenzOde, 0 10, 20 5 -5, options);01020304050-20-1001020-30-20-100102030MATLAB 程式設計進階篇：常微分方程式提示n要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入lorenz 指令 MATLAB 程式設計進階篇：常微分方程式範例-4 (I)n假設 OutputFcn 設成“myfunc”： n options

25、= odeset(OutputFcn, myfunc)nODE 指令會呼叫 myfunc(tspan, y0, init) 讓 myfunc 進行各種初始化動作 n積分步驟中，ODE 指令會持續呼叫status=myfunc(t, y)n假设 status=1，則停止積分 n積分結束時，ODE 指令會呼叫 myfunc( , , done)，讓 myfunc 進行收尾動作 nOutputSel 可指定要傳送到 OutputFcn 的狀態變數之元素 MATLAB 程式設計進階篇：常微分方程式範例-4 (II)n只要傳送第一及第三個 Lorenz 渾沌方程式的狀態變數至 odeplot n範例11

26、-9：odeOutputSelect01.m options = odeset(OutputSel, 1,3) % 只畫出第一和第三個狀態變數ode45(lorenzOde, 0 10, 20 5 -5, options);012345678910-30-20-1001020304050MATLAB 程式設計進階篇：常微分方程式範例-5 (I)nRefine 欄位可以使用內差法來增加輸出狀態變數的密度，以得到較平滑的輸出曲線 n用 Refine 欄位使 ode23 的輸出點個數增為原先三倍：n範例11-10：odeRefine01.m subplot(2,1,1);ode23(vdp1, 0

27、25, 3 3);subplot(2,1,2);options = odeset(Refine, 3);ode23(vdp1, 0 25, 3 3, options);MATLAB 程式設計進階篇：常微分方程式範例-5 (II)0510152025-4-20240510152025-4-2024MATLAB 程式設計進階篇：常微分方程式範例-6 n當 Stat=on 時，ODE 指令會在執行完畢後顯示計算過程的各種統計數字 n範例11-11：odeShowStats01.mn 71 successful stepsn 10 failed attemptsn 487 function evalu

28、ationst, y = ode45(vdp1, 0 25, 3 3, odeset(Stat, on);MATLAB 程式設計進階篇：常微分方程式範例-7n相同的統計數字，也可由 ODE 指令的第三個輸出變數傳回 n範例11-12：odeShowStats02.mn s =n 71n 10n 487n 0n 0n 0 t, y, s = ode45(vdp1, 0 25, 3 3);sMATLAB 程式設計進階篇：常微分方程式說明nMaxStep 及 InitialStep 欄位可用來調整最大積分步長及起始積分步長 n一般而言，不必去調整這兩個數值，因為 ODE 指令本身就具有步長自動調適功

29、能n若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加 n如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟 MATLAB 程式設計進階篇：常微分方程式11-4 ODE 檔案的進階用法n更進一步介紹 ODE 檔案的進階用法，使 ODE 指令能夠迅速且準確地算出積分結果 n可將 tspan積分時間範圍）、y0起始值及 optionsODE參數置於 ODE 檔案中，這些變數必須能由 ODE 檔案傳回，其格式為：n tspan, y0, option

30、s = odeFile(, , init) n假設 odeFile 即是我們的 ODE 檔案且 odeFile 滿足上述要求，則可以直接呼叫 ODE 指令如下：n t, y = solver(odeFile) n其中 solver 為前述的任一個 ODE 指令，它可由 odeFile 直接得到 tspan、y0 及 options 等積分所需的資訊 MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-1 (I)n以前述的 van der Pol 為例，若要能夠傳回 tspan、y0 及 options，vdp1.m 須改寫如下vdp3.m）：n type vdp3.m n f

31、unction output1, output2, output3 = vdp3(t, y, flag)n if strcmp(flag, )n mu = 1;n output1 = y(2); mu*(1-y(1)2)*y(2)-y(1);% dy/dtn elseif strcmp(flag, init),n output1 = 0; 25;% Time spann output2 = 3; 3;% Initial conditionsn output3 = odeset;% ODE optionsn end MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-1 (II)n範例11-13：odeAdvanced01.m ode45(vdp3)0510152025-3-2-101234MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-2 (I) nvan der Pol 的微分方程式有一個參數