常微分方程式课件

MATLAB程序设计进阶篇

常微分方程式張智星清大资工系多媒体检索实验室11-1ODE指令列表MATLAB用于求解常微分方程式的指令:

指令方法適用ODE類別ode45ExplicitRunge-Kutta(4,5)pairofDormand-PrinceNonstiffODEode23ExplicitRunge-Kutta(2,3)pairofBogackiandShampineNonstiffODEode113VariableorderAdams-Bashforth-MoultonPECEsolverNonstiffODEode15sNumericaldifferentiationformulas(NDFS)StiffODEode23sModifiedRosenbrockformulaoforder2StiffODEode23tTrapezoidalrulewitha“free”interpolantStiffODEode23tbImplicitRunge-KuttaformulawithabackwarddifferentiationformulaofordertwoStiffODEODE指令列表指令项目繁多，最主要可分两大类适用于Nonstiff系统一般的常微分方程式都是Nonstiff系统直接选用ode45、ode23或ode113来求解适用Stiff系统速率（即微分值）差异相常大使用一般的ode45、ode23或ode113来求解，可能会使得积分的步长（StepSizes）变得很小，以便降低积分误差至可容忍范围以内，会导致计算时间过长专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode23tb11-2ODE指令基本用法使用ODE指令时，必须先将要求解的ODE表示成一个函式输入为t（时间）及y（状态变量，StateVariables）输出则为dy（状态变量的微分值）ODE函式的档名为odeFile.m，则呼叫ODE指令的格式如下：[t,y]=solver('odeFile',[t0,t1],y0);[t0,t1]是积分的时间区间y0代表起始条件（InitialConditions）solver是前表所列的各种ODE指令t是输出的时间向量y则是对应的状态变量向量ODE指令基本用法以vanderPol微分方程为例，其方程式为：化成标准格式可用向量来表示成一般化的数学式为一向量，代表状态变量ODE指令基本用法假设=1，ODE档案（vdp1.m）可显示如下：>>typevdp1.m

functiondy=vdp1(t,y)mu=1;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];有了ODE档案，即可选用前述之ODE指令来求解ODE指令基本用法範例-1(II)

[0,25]代表积分的时间区间，[33]’

则代表起始条件（必须以行向量来表示）因为没有输出变量，所以上述程序执行结束后，MATLAB只会画出状态变量对时间的图形ODE指令基本用法範例-2(I)要取得积分所得的状态变量及对应的时间，可以加上输出变量以取得这些数据范例11-2：odeGetData01.m[t,y]=ode45('vdp1',[025],[33]');plot(t,y(:,1),t,y(:,2),':');xlabel('Timet');ylabel('Solutiony(t)andy''(t)');legend('y(t)','y''(t)');ODE指令基本用法範例-2(II)ODE指令基本用法範例-3(II)当值越来越大时，vanderPol方程式就渐渐变成一个Stiff的系统，此时若要解此方程式，就必须改用专门对付Stiff系统的指令ODE指令基本用法範例-4(I)将值改成1000，ODE档案改写成（vdp2.m）：

>>typevdp2.mfunctiondy=vdp2(t,y)mu=1000;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];ODE指令基本用法範例-4(II)选用专门对付Stiff系统的ODE指令，例如ode15s，来求解此系统并作图显示範例11-4：ode15s01.m[t,y]=ode15s('vdp2',[03000],[21]');subplot(2,1,1);plot(t,y(:,1),'-o');xlabel('Timet');ylabel('y(t)');subplot(2,1,2);plot(t,y(:,2),'-o');xlabel('Timet');ylabel('y''(t)'); %注意單引號「'」的重覆使用ODE指令基本用法範例-5(I)若要画出二维平面相位图，可以使用下列范例：范例11-5：ode15s02.m若要产生在某些特定时间点的状态变量值，则呼叫ODE指令的格式可改成：[t,y]=solver('odeFile',[t0,t1,…,tn],y0);其中[t0,t1,…,tn]即是特定时间点所形成的向量[t,y]=ode15s('vdp2',[03000],[21]');subplot(1,1,1);plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)')ODE指令基本用法範例-5(II)11-3ODE指令的選項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指令的選項>>odesetAbsTol:[positivescalarorvector{1e-6}]RelTol:[positivescalar{1e-3}]NormControl:[on|{off}]NonNegative:[vectorofintegers]OutputFcn:[function_handle]OutputSel:[vectorofintegers]Refine:[positiveinteger]Stats:[on|{off}]InitialStep:[positivescalar]MaxStep:[positivescalar]BDF:[on|{off}]MaxOrder:[1|2|3|4|{5}]Jacobian:[matrix|function_handle]JPattern:[sparsematrix]Vectorized:[on|{off}]Mass:[matrix|function_handle]MStateDependence:[none|{weak}|strong]MvPattern:[sparsematrix]MassSingular:[yes|no|{maybe}]InitialSlope:[vector]Events:[function_handle]由odeset產生的ODE選項類別欄位名稱資料型態預設值說明誤差容忍度之相關欄位RelTol正純量相對誤差容忍度AbsTo1正純量或向量絕對誤差容忍度積分輸出之相關欄性OutPutFcn字串‘odeplot’輸出函式（若ODE指令無輸出變數，則在數值積分執行完畢後，MATLAB會呼叫此輸出函式）OutputSel索引向量全部ODE指令之輸出變數的索引值，以決定那些輸出變數之元素將被送到輸出函式Refine正整數1或4(forode45)Refine=2可產生兩倍數量的輸出點，Refine=3可產生三倍數量的輸出點，依此類推。Statson或offoffStats='on'產生計算過程的各種統計資料。由odeset產生的ODE選項若F(t,y,’Jacobian')傳回若F([],[],’JPattern’)傳回，且Jacobian矩陣之相關欄位Jconstanton或offoff如果Jacobian矩陣常數，則JConstant='on'Jacobianon或offoff，則Jacobian='on‘Jpatternon或offoff是稀疏矩陣，則JPattem='on'Vectorizedon或offoff若F(t,[y1,y2…..])傳回[F(t,y1),F(t,y2)…..]，則Vectorized='on'積分步長（StepSizes）之相關欄位MaxStep正純量ODE指令之積分步長的上限InitialStep正純量起始步長的建議值常用到的欄位來進行說明f在积分误差容忍度方面，每一次积分所产生的局部误差e(i)，必须满足下列方程式：max(RelTol*,AbsTol(i))其中i代表第i个状态变量降低RelTol及AbsTol来求得更精确的积分结果範例-1(I)範例11-6：odeRelTol01.msubplot(2,1,1);ode45('vdp1',[025],[33]');title('RelTol=0.01');options=odeset('RelTol',0.00001);subplot(2,1,2);ode45('vdp1',[025],[33]',options);title('RelTol=0.0001');積分輸出方面說明以Lorenz渾沌方程式（LorenzChaoticEquation）為例

>>typelorenzOde.mfunctiondy=lorenzOde(t,y)%LORENZODE:ODELorenzchaoticequationIGMA=10.;RHO=28;BETA=8./3.;A=[-BETA0y(2)0-SIGMASIGMA-y(2)RHO-1];dy=A*y;範例-2使用

ode45對Lorenz渾沌方程式進行數值積分範例11-7：odeLorenz01.m上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形ode45('lorenzOde',[010],[205-5]');範例-3上述範例畫三度空間之相位圖形範例11-8：odeLorenz02.m

圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線options=odeset('OutputFcn','odephas3');%使用odephas3進行繪圖ode45('lorenzOde',[010],[205-5]',options);提示要觀看Lorenz渾沌方程式隨時間而變的動畫，可在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的狀態變數之元素範例-4(II)只要傳送第一及第三個Lorenz渾沌方程式的狀態變數至odeplot範例11-9：odeOutputSelect01.moptions=odeset('OutputSel',[1,3])%只畫出第一和第三個狀態變數ode45('lorenzOde',[010],[205-5]',options);範例-5(I)Refine欄位可以使用內差法來增加輸出狀態變數的密度，以得到較平滑的輸出曲線用Refine欄位使ode23的輸出點個數增為原先三倍：範例11-10：odeRefine01.msubplot(2,1,1);ode23('vdp1',[025],[33]');subplot(2,1,2);options=odeset('Refine',3);ode23('vdp1',[025],[33]',options);範例-5(II)範例-6當Stat=on時，ODE指令會在執行完畢後顯示計算過程的各種統計數字範例11-11：odeShowStats01.m71successfulsteps10failedattempts487functionevaluations[t,y]=ode45('vdp1',[025],[33]',odeset('Stat','on'));範例-7相同的統計數字，也可由ODE指令的第三個輸出變數傳回

範例11-12：odeShowStats02.ms=7110487000[t,y,s]=ode45('vdp1',[025],[33]');s說明MaxStep及InitialStep欄位可用來調整最大積分步長及起始積分步長一般而言，不必去調整這兩個數值，因為ODE指令本身就具有步長自動調適功能若要產生更多輸出點，可直接調整Refine欄位值。調整MaxStep雖然可以達到同樣效果，但是計算時間可能會大幅增加如果積分結果不甚準確，請勿先調降MaxStep，您應先調降RelTol及AbsTol。調降MaxStep是最後的步驟11-4ODE檔案的進階用法更進一步介紹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等積分所需的資訊ODE檔案的進階用法範例-1(I)以前述的vanderPol為例，若要能夠傳回tspan、y0及options，vdp1.m須改寫如下（vdp3.m）：>>typevdp3.mfunction[output1,output2,output3]=vdp3(t,y,flag)ifstrcmp(flag,'') mu=1; output1=[y(2);mu*(1-y(1)^2)*y(2)-y(1)]; %dy/dtelseifstrcmp(flag,'init'), output1=[0;25]; %Timespan output2=[3;3]; %Initialconditions output3=odeset; %ODEoptionsendODE檔案的進階用法範例-1(II)範例11-13：odeAdvanced01.mode45('vdp3')ODE檔案的進階用法範例-2(I)vanderPol的微分方程式有一個參數，希望從外面傳入此參數的值(vdp4.m)>>typevdp4.mfunction[output1,output2,output3]=vdp4(t,y,flag,mu)ifnargin<4|isempty(mu), mu=1;endifstrcmp(f