MATLAB程序設(shè)計(jì)進(jìn)階篇

常微分方程式張智星清大資工系多媒體檢索實(shí)驗(yàn)室11-1ODE指令列表MATLAB用于求解常微分方程式的指令:

指令方法適用ODE類別ode45ExplicitRunge-Kutta(4,5)pairofDormand-PrinceNonstiffODEode23ExplicitRunge-Kutta(2,3)pairofBogackiandShampineNonstiffODEode113VariableorderAdams-Bashforth-MoultonPECEsolverNonstiffODEode15sNumericaldifferentiationformulas(NDFS)StiffODEode23sModifiedRosenbrockformulaoforder2StiffODEode23tTrapezoidalrulewitha“free”interpolantStiffODEode23tbImplicitRunge-KuttaformulawithabackwarddifferentiationformulaofordertwoStiffODEODE指令列表指令項(xiàng)目繁多，最主要可分兩大類適用于Nonstiff系統(tǒng)一般的常微分方程式都是Nonstiff系統(tǒng)直接選用ode45、ode23或ode113來求解適用Stiff系統(tǒng)速率（即微分值）差異相常大使用一般的ode45、ode23或ode113來求解，可能會(huì)使得積分的步長（StepSizes）變得很小，以便降低積分誤差至可容忍范圍以內(nèi)，會(huì)導(dǎo)致計(jì)算時(shí)間過長專門對付Stiff系統(tǒng)的指令，例如ode15s、ode23s、ode23t及ode23tb11-2ODE指令基本用法使用ODE指令時(shí)，必須先將要求解的ODE表示成一個(gè)函式輸入為t（時(shí)間）及y（狀態(tài)變量，StateVariables）輸出則為dy（狀態(tài)變量的微分值）ODE函式的檔名為odeFile.m，則呼叫ODE指令的格式如下：[t,y]=solver('odeFile',[t0,t1],y0);[t0,t1]是積分的時(shí)間區(qū)間y0代表起始條件（InitialConditions）solver是前表所列的各種ODE指令t是輸出的時(shí)間向量y則是對應(yīng)的狀態(tài)變量向量ODE指令基本用法以vanderPol微分方程為例，其方程式為：化成標(biāo)準(zhǔn)格式可用向量來表示成一般化的數(shù)學(xué)式為一向量，代表狀態(tài)變量ODE指令基本用法假設(shè)=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指令基本用法範(fàn)例-1(II)

[0,25]代表積分的時(shí)間區(qū)間，[33]’

則代表起始條件（必須以行向量來表示）因?yàn)闆]有輸出變量，所以上述程序執(zhí)行結(jié)束后，MATLAB只會(huì)畫出狀態(tài)變量對時(shí)間的圖形ODE指令基本用法範(fàn)例-2(I)要取得積分所得的狀態(tài)變量及對應(yīng)的時(shí)間，可以加上輸出變量以取得這些數(shù)據(jù)范例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指令基本用法範(fàn)例-2(II)ODE指令基本用法範(fàn)例-3(II)當(dāng)值越來越大時(shí)，vanderPol方程式就漸漸變成一個(gè)Stiff的系統(tǒng)，此時(shí)若要解此方程式，就必須改用專門對付Stiff系統(tǒng)的指令ODE指令基本用法範(fàn)例-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指令基本用法範(fàn)例-4(II)選用專門對付Stiff系統(tǒng)的ODE指令，例如ode15s，來求解此系統(tǒng)并作圖顯示範(fàn)例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)'); %注意單引號(hào)「'」的重覆使用ODE指令基本用法範(fàn)例-5(I)若要畫出二維平面相位圖，可以使用下列范例：范例11-5：ode15s02.m若要產(chǎn)生在某些特定時(shí)間點(diǎn)的狀態(tài)變量值，則呼叫ODE指令的格式可改成：[t,y]=solver('odeFile',[t0,t1,…,tn],y0);其中[t0,t1,…,tn]即是特定時(shí)間點(diǎn)所形成的向量[t,y]=ode15s('vdp2',[03000],[21]');subplot(1,1,1);plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)')ODE指令基本用法範(fàn)例-5(II)11-3ODE指令的選項(xiàng)ODE指令可以接受第四個(gè)輸入變量，代表積分過程用到的各種選項(xiàng)（Options），此種ODE指令的格式為：[t,y]=solver('odeFile',[t0,tn],y0,options);其中options是由odeset指令來控制的結(jié)構(gòu)變量結(jié)構(gòu)變量即包含了積分過程用到的各種選項(xiàng)odeset的一般格式如下：options=odeset('name1',value1,'name2',value2,…)其字段name1的值為value1，字段name2的值為value2，依此類推未被設(shè)定的字段，其域值即為默認(rèn)值ODE指令的選項(xiàng)>>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產(chǎn)生的ODE選項(xiàng)類別欄位名稱資料型態(tài)預(yù)設(shè)值說明誤差容忍度之相關(guān)欄位RelTol正純量相對誤差容忍度AbsTo1正純量或向量絕對誤差容忍度積分輸出之相關(guān)欄性O(shè)utPutFcn字串‘odeplot’輸出函式（若ODE指令無輸出變數(shù)，則在數(shù)值積分執(zhí)行完畢後，MATLAB會(huì)呼叫此輸出函式）OutputSel索引向量全部ODE指令之輸出變數(shù)的索引值，以決定那些輸出變數(shù)之元素將被送到輸出函式Refine正整數(shù)1或4(forode45)Refine=2可產(chǎn)生兩倍數(shù)量的輸出點(diǎn)，Refine=3可產(chǎn)生三倍數(shù)量的輸出點(diǎn)，依此類推。Statson或offoffStats='on'產(chǎn)生計(jì)算過程的各種統(tǒng)計(jì)資料。由odeset產(chǎn)生的ODE選項(xiàng)若F(t,y,’Jacobian')傳回若F([],[],’JPattern’)傳回，且Jacobian矩陣之相關(guān)欄位Jconstanton或offoff如果Jacobian矩陣常數(shù)，則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）之相關(guān)欄位MaxStep正純量ODE指令之積分步長的上限InitialStep正純量起始步長的建議值常用到的欄位來進(jìn)行說明f在積分誤差容忍度方面，每一次積分所產(chǎn)生的局部誤差e(i)，必須滿足下列方程式：max(RelTol*,AbsTol(i))其中i代表第i個(gè)狀態(tài)變量降低RelTol及AbsTol來求得更精確的積分結(jié)果範(fàn)例-1(I)範(fàn)例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;範(fàn)例-2使用

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

圖形中只出現(xiàn)一條曲線，此曲線代表以三個(gè)狀態(tài)變數(shù)為座標(biāo)、以時(shí)間為參數(shù)的一條三度空間中的曲線options=odeset('OutputFcn','odephas3');%使用odephas3進(jìn)行繪圖ode45('lorenzOde',[010],[205-5]',options);提示要觀看Lorenz渾沌方程式隨時(shí)間而變的動(dòng)畫，可在MATLAB指令視窗下直接輸入lorenz指令範(fàn)例-4(I)假設(shè)OutputFcn設(shè)成“myfunc”：

options=odeset('OutputFcn','myfunc')ODE指令會(huì)呼叫myfunc(tspan,y0,‘init’)讓myfunc進(jìn)行各種初始化動(dòng)作積分步驟中，ODE指令會(huì)持續(xù)呼叫status=myfunc(t,y)若status=1，則停止積分積分結(jié)束時(shí)，ODE指令會(huì)呼叫myfunc([],[],‘done’)，讓myfunc進(jìn)行收尾動(dòng)作OutputSel可指定要傳送到OutputFcn的狀態(tài)變數(shù)之元素範(fàn)例-4(II)只要傳送第一及第三個(gè)Lorenz渾沌方程式的狀態(tài)變數(shù)至odeplot範(fàn)例11-9：odeOutputSelect01.moptions=odeset('OutputSel',[1,3])%只畫出第一和第三個(gè)狀態(tài)變數(shù)ode45('lorenzOde',[010],[205-5]',options);範(fàn)例-5(I)Refine欄位可以使用內(nèi)差法來增加輸出狀態(tài)變數(shù)的密度，以得到較平滑的輸出曲線用Refine欄位使ode23的輸出點(diǎn)個(gè)數(shù)增為原先三倍：範(fàn)例11-10：odeRefine01.msubplot(2,1,1);ode23('vdp1',[025],[33]');subplot(2,1,2);options=odeset('Refine',3);ode23('vdp1',[025],[33]',options);範(fàn)例-5(II)範(fàn)例-6當(dāng)Stat=on時(shí)，ODE指令會(huì)在執(zhí)行完畢後顯示計(jì)算過程的各種統(tǒng)計(jì)數(shù)字範(fàn)例11-11：odeShowStats01.m71successfulsteps10failedattempts487functionevaluations[t,y]=ode45('vdp1',[025],[33]',odeset('Stat','on'));範(fàn)例-7相同的統(tǒng)計(jì)數(shù)字，也可由ODE指令的第三個(gè)輸出變數(shù)傳回

範(fàn)例11-12：odeShowStats02.ms=7110487000[t,y,s]=ode45('vdp1',[025],[33]');s說明MaxStep及InitialStep欄位可用來調(diào)整最大積分步長及起始積分步長一般而言，不必去調(diào)整這兩個(gè)數(shù)值，因?yàn)镺DE指令本身就具有步長自動(dòng)調(diào)適功能若要產(chǎn)生更多輸出點(diǎn)，可直接調(diào)整Refine欄位值。調(diào)整MaxStep雖然可以達(dá)到同樣效果，但是計(jì)算時(shí)間可能會(huì)大幅增加如果積分結(jié)果不甚準(zhǔn)確，請勿先調(diào)降MaxStep，您應(yīng)先調(diào)降RelTol及AbsTol。調(diào)降MaxStep是最後的步驟11-4ODE檔案的進(jìn)階用法更進(jìn)一步介紹ODE檔案的進(jìn)階用法，使ODE指令能夠迅速且準(zhǔn)確地算出積分結(jié)果可將tspan（積分時(shí)間範(fàn)圍）、y0（起始值）及options（ODE參數(shù)）置於ODE檔案中，這些變數(shù)必須能由ODE檔案傳回，其格式為：[tspan,y0,options]=odeFile([],[],'init')假設(shè)odeFile即是我們的ODE檔案且odeFile滿足上述要求，則可以直接呼叫ODE指令如下：[t,y]=solver('odeFile')其中solver為前述的任一個(gè)ODE指令，它可由odeFile直接得到tspan、y0及options等積分所需的資訊ODE檔案的進(jìn)階用法範(fàn)例-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檔案的進(jìn)階用法範(fàn)例-1(II)範(fàn)例11-13：odeAdvanced01.mode45('vdp3')ODE檔案的進(jìn)階用法範(fàn)例-2(I)vanderPol的微分方程式有一個(gè)參數(shù)，希望從外面?zhèn)魅氪藚?shù)的值(vdp4.m)>>typevdp4.mfunction[output1,output2,output3]=vdp4(t,y,flag,mu)ifnargin<4|isempty(mu), mu=1;endifstrcmp(f