《《常微分方程式》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《常微分方程式》PPT课件.ppt（46页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、MATLAB 程式設計進階篇：常微分方程式MATLAB 程序设计进阶篇常微分方程式張智星http:/www.cs.nthu.edu.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 ODEode113Variab
2、le order 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 o
3、rder 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 程式設計進階篇：常微分方程式提示
4、n使用 Simulink 來求解常微分方程式nSimulink是和MATLAB共同使用的一套软件n可使用拖拉的方式来建立动态系统n可直接产生C程序代码或进行动画显示n功能非常强大 MATLAB 程式設計進階篇：常微分方程式11-2 ODE 指令基本用法 n使用ODE指令时，必须先将要求解的ODE表示成一个函式n输入为t（时间）及y（状态变量，State Variables）n输出则为dy（状态变量的微分值）nODE函式的档名为odeFile.m，则呼叫ODE指令的格式如下：t,y=solver(odeFile,t0,t1,y0);nt0,t1是积分的时间区间ny0代表起始条件（Initial
5、Conditions）nsolver是前表所列的各种ODE指令nt是输出的时间向量ny则是对应的状态变量向量MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法n以van der Pol微分方程为例，其方程式为：n化成标准格式n可用向量来表示成一般化的数学式n 为一向量，代表状态变量MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法n假设=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);n有了ODE档案，即可选用前述之ODE指令来求解MATL
6、AB 程式設計進階篇：常微分方程式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只会画出状态变量对时间的图形MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-2(I)n要取得积分所得的状态变量及对应的
7、时间，可以加上输出变量以取得这些数据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)MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-3(I)n也可以画出 及 在 相位平面（Phase Plane）的运动情况n範例11-3：odePhastPlot01.mt,y=ode45(vdp1
8、,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的系统，此时若要解此方程式，就必须改用专门对付Stiff系统的指令MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-4(I)n将 值改成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
9、(1);MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-4(II)n选用专门对付Stiff系统的ODE指令，例如ode15s，来求解此系统并作图显示n範例11-4：ode15s01.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
10、)n由上图可知，的变化相当剧烈（超过 ），这就是Stiff系统的特色MATLAB 程式設計進階篇：常微分方程式ODE 指令基本用法範例-5(I)n若要画出二维平面相位图，可以使用下列范例：n范例11-5：ode15s02.mn若要产生在某些特定时间点的状态变量值，则呼叫ODE指令的格式可改成：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),y(:,2),-o);xlabel(y(t);ylabel(y(t)MATLAB 程
11、式設計進階篇：常微分方程式ODE 指令基本用法範例-5(II)MATLAB 程式設計進階篇：常微分方程式11-3 ODE 指令的選項nODE指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为：t,y=solver(odeFile,t0,tn,y0,options);n其中options是由odeset指令来控制的结构变量n结构变量即包含了积分过程用到的各种选项nodeset的一般格式如下：options=odeset(name1,value1,name2,value2,)n其字段name1的值为value1，字段name2的值为value2，依此类
12、推n未被设定的字段，其域值即为默认值MATLAB 程式設計進階篇：常微分方程式ODE 指令的選項n也可以只改变一个现存的options结构变量中，某个字段的值，其格式如下：newOptions=odeset(options,name,value);n若要读出某个字段的值，可用odeget，其格式如下：value=odeget(otpions,name);n其中name为域名，value即为对应的域值n当使用odeset指令时，若无任何输入变量，则odeset会显示所有的域名及域值，并以大括号代表默认值MATLAB 程式設計進階篇：常微分方程式ODE 指令的選項 odeset AbsTol:po
13、sitive 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|fun
14、ction_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 MATLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項 類別欄位名稱資料型態預設值說明誤差容忍度之相關欄位RelTol正純量相對誤差容忍度AbsTo1正純量
15、或向量絕對誤差容忍度積分輸出之相關欄性OutPutFcn字串odeplot輸出函式（若 ODE 指令無輸出變數，則在數值積分執行完畢後，MATLAB 會呼叫此輸出函式）OutputSel索引向量全部ODE 指令之輸出變數的索引值，以決定那些輸出變數之元素將被送到輸出函式Refine正整數1或4(for ode45)Refine=2 可產生兩倍數量的輸出點，Refine=3 可產生三倍數量的輸出點，依此類推。Statson 或 offoffStats=on 產生計算過程的各種統計資料。MATLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項若F(t,y,Jacobian)
16、傳回 若 F(,JPattern)傳回，且Jacobian 矩陣之相關欄位Jconstanton 或 offoff如果 Jacobian 矩陣常數，則 JConstant=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正純量起始步長的建議值MA
17、TLAB 程式設計進階篇：常微分方程式由 odeset 產生的 ODE 選項質量矩陣之相關欄位Massnone，M，M(t)，或 M(t,y)none表明 ODE 指令案是否會傳回質量矩陣MassSingular yes，no 或 maybemaybe表明質量矩陣是否為 Singular事件發生時間之相關欄位Eventson 或 offoff若 ODE 檔案並傳回事件函式及相關資訊，則 Events=onOde15s 之相關欄位MaxOrder1,2,3,4 或 55積分公式用到的最高次數BDFon 或 offoff若使用 BDF（Backward Differentiation Formul
18、a）則 BDF=onMATLAB 程式設計進階篇：常微分方程式常用到的欄位來進行說明 nf在积分误差容忍度方面，每一次积分所产生的局部误差e(i)，必须满足下列方程式：max(RelTol*,AbsTol(i)n其中i代表第i个状态变量n降低RelTol及AbsTol来求得更精确的积分结果MATLAB 程式設計進階篇：常微分方程式範例-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
19、(vdp1,0 25,3 3,options);title(RelTol=0.0001);MATLAB 程式設計進階篇：常微分方程式範例-1(II)n第一個圖所使用的相對誤差值是0.01（預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長 MATLAB 程式設計進階篇：常微分方程式積分輸出方面說明n积分输出的相关处理方面n选用一个OutputFcnn当ODE指令没有输出变量时，此输出函式OutputFcn会被MATLAB呼叫nOutputFcn的默认值是”odeplot”，其功能为画出所有的状态变量n其它可用的函式nodephas2：画出2
20、-D的相位平面（Phase Plane）nodephas3：画出3-D的相位平面nodeprint：随时在指令窗口印出计算结果MATLAB 程式設計進階篇：常微分方程式積分輸出方面說明n以 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*
21、y;MATLAB 程式設計進階篇：常微分方程式範例-2 n使用 ode45 對Lorenz 渾沌方程式進行數值積分n範例11-7：odeLorenz01.mn上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形 ode45(lorenzOde,0 10,20 5-5);MATLAB 程式設計進階篇：常微分方程式範例-3n上述範例畫三度空間之相位圖形 n範例11-8：odeLorenz02.m n圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線 options=odeset(OutputFcn,odephas3);%使用 odephas3 進行繪圖od
22、e45(lorenzOde,0 10,20 5-5,options);MATLAB 程式設計進階篇：常微分方程式提示n要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入lorenz 指令 MATLAB 程式設計進階篇：常微分方程式範例-4(I)n假設 OutputFcn 設成“myfunc”：options=odeset(OutputFcn,myfunc)nODE 指令會呼叫 myfunc(tspan,y0,init)讓 myfunc 進行各種初始化動作 n積分步驟中，ODE 指令會持續呼叫status=myfunc(t,y)n若 status=1，則停止
23、積分 n積分結束時，ODE 指令會呼叫 myfunc(,done)，讓 myfunc 進行收尾動作 nOutputSel 可指定要傳送到 OutputFcn 的狀態變數之元素 MATLAB 程式設計進階篇：常微分方程式範例-4(II)n只要傳送第一及第三個 Lorenz 渾沌方程式的狀態變數至 odeplot n範例11-9：odeOutputSelect01.m options=odeset(OutputSel,1,3)%只畫出第一和第三個狀態變數ode45(lorenzOde,0 10,20 5-5,options);MATLAB 程式設計進階篇：常微分方程式範例-5(I)nRefine
24、欄位可以使用內差法來增加輸出狀態變數的密度，以得到較平滑的輸出曲線 n用 Refine 欄位使 ode23 的輸出點個數增為原先三倍：n範例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);MATLAB 程式設計進階篇：常微分方程式範例-5(II)MATLAB 程式設計進階篇：常微分方程式範例-6 n當 Stat=on 時，ODE 指令會在執行完畢後顯示計算過程的各種統計數字 n範例11-11：o
25、deShowStats01.m 71 successful steps 10 failed attempts 487 function evaluationst,y=ode45(vdp1,0 25,3 3,odeset(Stat,on);MATLAB 程式設計進階篇：常微分方程式範例-7n相同的統計數字，也可由 ODE 指令的第三個輸出變數傳回 n範例11-12：odeShowStats02.m s=71 10 487 0 0 0 t,y,s=ode45(vdp1,0 25,3 3);sMATLAB 程式設計進階篇：常微分方程式說明nMaxStep 及 InitialStep 欄位可用來調整最
26、大積分步長及起始積分步長 n一般而言，不必去調整這兩個數值，因為 ODE 指令本身就具有步長自動調適功能n若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加 n如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟 MATLAB 程式設計進階篇：常微分方程式11-4 ODE 檔案的進階用法n更進一步介紹 ODE 檔案的進階用法，使 ODE 指令能夠迅速且準確地算出積分結果 n可將 tspan（積分時間範圍）、y0（起始值）及 options（ODE
27、參數）置於 ODE 檔案中，這些變數必須能由 ODE 檔案傳回，其格式為：tspan,y0,options=odeFile(,init)n假設 odeFile 即是我們的 ODE 檔案且 odeFile 滿足上述要求，則可以直接呼叫 ODE 指令如下：t,y=solver(odeFile)n其中 solver 為前述的任一個 ODE 指令，它可由 odeFile 直接得到 tspan、y0 及 options 等積分所需的資訊 MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-1(I)n以前述的 van der Pol 為例，若要能夠傳回 tspan、y0 及 option
28、s，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 MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-1(II)n
29、範例11-13：odeAdvanced01.m ode45(vdp3)MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-2(I)nvan 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),o
30、utput1=0;25;%Time span output2=3;3;%Initial conditions output3=odeset;%ODE options end MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-2(II)n 就變成一個選擇性（Optional）的參數，其預設值為 1 n將 的值從 MATLAB 傳入，並畫出不同 值下的 van der Pol 方程式的狀態變數：n範例11-14：odeAdvanced02.msubplot(2,1,1);ode45(vdp4,1);%mu=1subplot(2,1,2);ode45(vdp4,3);%mu=3M
31、ATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法範例-2(III)n 的值分別是 1 及 3 n用到了許多空矩陣，這些空矩陣代表取用預設值，因此 ode45 會直接從 vdp4.m 取用時間區間及變數起始值 n也可以傳入二個或更多的參數，MATLAB 及 ODE 指令對於可傳入的參數個數並無設限 MATLAB 程式設計進階篇：常微分方程式ODE 檔案的進階用法列表n為解決其它較複雜的 ODEs 及 DAEs（Differential Algebra Equations），ODE 檔案亦可在不同的旗號（Flag）下傳回不同的資訊，列表如下：旗號傳回值（空字串）dy(=F(t,y)inittspan,yo 及 optionsjacobianJacobian 矩陣 J(t,y)=2F/2YjpatternJacobian sparsity pattern 之矩陣mass解 M(t,y)y=F(t,y)所須的 質量矩陣 Mevents定義事件發生點的各種資訊