Matlab在常微分方程求解中應(yīng)用

Matlab在常微分方程求解中應(yīng)用第一頁，共51頁。實(shí)驗(yàn)?zāi)康?/p>

（1）學(xué)會用Matlab軟件求解微分方程的初值問題（2）了解微分方程數(shù)值解思想，掌握基本的微分方程數(shù)值解方法（3）學(xué)會根據(jù)實(shí)際問題建立簡單微分方程數(shù)學(xué)模型（4）了解計(jì)算機(jī)數(shù)據(jù)仿真、數(shù)據(jù)模擬的基本方法

第二頁，共51頁。17世紀(jì)：初等解法18世紀(jì)：初等解法和無窮級數(shù)方法19世紀(jì)：解的存在性、奇點(diǎn)理論、定性理論、穩(wěn)定性理論

包含一個自變量和它的未知函數(shù)以及未知函數(shù)的導(dǎo)數(shù)的等式

形成和發(fā)展與力學(xué)、天文學(xué)、物理學(xué)及其他自然科學(xué)技術(shù)的發(fā)展互相促進(jìn)和推動

常微分方程第三頁，共51頁。定理

設(shè)函數(shù)在區(qū)域上連續(xù)，且在區(qū)域D內(nèi)滿足李普希茲(Lipschitz)條件，即存在正數(shù)L，使得對于R內(nèi)任意兩點(diǎn)與,恒有則初值問題(1)的解存在并且唯一。第四頁，共51頁。常微分方程的解析解

求微分方程（組）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)

結(jié)果：u=tg(t-c)記號：在表達(dá)微分方程時，用字母D表示求微分，D2、D3等表示求高階微分。D后所跟字母為因變量，自變量可以指定或由系統(tǒng)規(guī)則選定為缺省。例如：微分方程可以表示為D2y=0.第五頁，共51頁。

解輸入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')結(jié)果為:y=3e-2xsin（5x）第六頁，共51頁。解輸入命令：

[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t')；x=simple(x)%將x化簡y=simple(y)z=simple(z)結(jié)果為：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=(c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

第七頁，共51頁。雖然說解析解是最精確的，但是實(shí)際問題中常要求研究常微分方程的數(shù)值解第八頁，共51頁。常微分方程的數(shù)值解

常微分方程中只有一些典型方程能求出初等解(用初等函數(shù)表示的解)。另外，有些初值問題雖然有初等解，但由于形式太復(fù)雜不便于應(yīng)用。因此，有必要探討常微分方程初值問題的數(shù)值解法。以下主要介紹一階常微分方程初值問題幾種經(jīng)典數(shù)值解方法：歐拉法及改進(jìn)的歐拉法。其它方法：龍格-庫塔法、阿達(dá)姆斯方法；一階微分方程組與高階方程初值問題的數(shù)值解法；二階常微分方程值問題的差分方法等。第九頁，共51頁。求解常微分方程初值問題的數(shù)值解的整體思路：

（1）尋求準(zhǔn)確解在一系列離散節(jié)點(diǎn)：上的近似值

稱為問題的數(shù)值解，數(shù)值解所滿足的離散方程統(tǒng)稱為差稱為步長，實(shí)用中常取定步長。分格式，第十頁，共51頁。建立數(shù)值解法的一些途徑1、用差商代替導(dǎo)數(shù)

若步長h較小，則有故有公式：此即歐拉法第十一頁，共51頁。Euler法的幾何意義:找到了積分的一條近似折線第十二頁，共51頁。2、使用數(shù)值積分對方程y`=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實(shí)際應(yīng)用時，與歐拉公式結(jié)合使用：此即改進(jìn)的歐拉法故有公式：第十三頁，共51頁。第十四頁，共51頁。第十五頁，共51頁。第十六頁，共51頁。3、使用泰勒公式

以此方法為基礎(chǔ)，有龍格-庫塔法、線性多步法等方法。4、數(shù)值公式的精度

當(dāng)一個數(shù)值公式的截?cái)嗾`差可表示為O（hk+1）時（k為正整數(shù)，h為步長），稱它是一個k階公式。k越大，則數(shù)值公式的精度越高。歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。第十七頁，共51頁。ODE指令列表

MATLAB用於求解常微分方程式的指令:

指令方法應(yīng)用ODE類別ode45ExplicitRunge-Kutta(4,5)pairofDormand-PrinceNonstiffODEode23ExplicitRunge-Kutta(2,3)pairofBogackiandShampineNonstiffODEode113VariableorderAdams-Bashforth-MoultonPECEsolverNonstiffODEode15sNumericaldifferentiationformulas(NDFS)StiffODEode23sModifiedRosenbrockformulaoforder2StiffODEode23tTrapezoidalrulewitha“free”interpolantStiffODEode23tbImplicitRunge-KuttaformulawithabackwarddifferentiationformulaofordertwoStiffODE第十八頁，共51頁。

適用于

Nonstiff系統(tǒng)

速率（即微分值）差異相常大使用一般的ode45、ode23或ode113來求解,可能會使得積分的步長(StepSizes)變得很小,以便降低積分誤差至可容忍范圍以內(nèi),會導(dǎo)致計(jì)算時間過長專門對付Stiff系統(tǒng)的指令，例如ode15s、ode23s、

ode23t及ode23tb指令項(xiàng)目繁多，最主要可分兩大類：適用于

Nonstiff系統(tǒng)

一般的常微分方程式都是Nonstiff系統(tǒng)直接采用ode45、ode23或ode113來求解第十九頁，共51頁。提示使用Simulink來求解常微分方程式Simulink是和MATLAB共同使用的一套軟件可使用拖拉的方式來建立動力系統(tǒng)可直接產(chǎn)生C語言代碼或進(jìn)行動畫演示功能非常強(qiáng)大第二十頁，共51頁。ODE指令基本用法

使用ODE指令時，必須先將要求解的ODE表示

成一個函數(shù)

輸入為t(時間)及y(狀態(tài)函數(shù):StateVariables）輸出則為dy（狀態(tài)變量的微分值）第二十一頁，共51頁。ODE指令基本用法

若ODE函數(shù)的文件為odeFile.m,則調(diào)用ODE指令

的格式如下：

[t,y]=solver('odeFile',[t0,t1],y0)

[t0,t1]是積分的時間區(qū)間y0代表起始條件（InitialConditions）solver是前表所列的各種ODE指令t是輸出的時間向量y是對應(yīng)的狀態(tài)變量向量第二十二頁，共51頁。[t,y]=solver（’f’,ts,y0,options）ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0,tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫塔-芬爾格算法ode45：運(yùn)用組合的4/5階龍格-庫塔-芬爾格算法自變量值函數(shù)值用于設(shè)定誤差限(缺省時設(shè)定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對誤差和絕對誤差.第二十三頁，共51頁。

1、在解n個未知函數(shù)的方程組時，y0和y均為n維向量，m-文件中的待解方程組應(yīng)以y的分量形式寫成.

2、使用Matlab軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.注意:第二十四頁，共51頁。VanderPol微分方程1928年荷蘭的范德波耳(VanderPol)為描述LC回路的電子管振蕩器建立了著名的vanderPol方程.它在自激振蕩理論中有著重要的意義,一直作為數(shù)學(xué)物理方程中的一個基本方程.這是一個具有可變非線性阻尼的微分方程,代表了一類極為典型的非線性問題.和其他非線性微分方程在數(shù)學(xué)上無法精確求解一樣,人們一直在努力尋找求解這類方程近似解析解的方法,并樂于用VanderPol方程來檢驗(yàn)求解方法的有效性.

第二十五頁，共51頁。VanderPol微分方程其方程形式為：

令

將VanderPol微分方程化成標(biāo)準(zhǔn)形式第二十六頁，共51頁。寫成向量的形式：

為一個向量，代表狀態(tài)變量第二十七頁，共51頁。設(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)];建立了vdp1.m后，即可選用前述ODE指令來求解第二十八頁，共51頁。在=1時，vanderPol方程并非Stiff系統(tǒng)，所以使用ode45來顯示結(jié)果

odeBasic01.mode45('vdp1',[025],[33]');

[025]代表積分的時間區(qū)間，[33]’則代表起始條件因?yàn)闆]有輸出變量，所以上述程序執(zhí)行結(jié)束后，MATLAB只畫出狀態(tài)變量對時間的圖形第二十九頁，共51頁。第三十頁，共51頁。

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)');第三十一頁，共51頁。第三十二頁，共51頁。也可以畫出及在相位平面（PhasePlane）的運(yùn)動情況

odePhastPlot01.m[t,y]=ode45('vdp1',[025],[33]');plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)');第三十三頁，共51頁。第三十四頁，共51頁。當(dāng)值越來越大時，VanderPol方程就漸變?yōu)橐粋€Stiff系統(tǒng)，此時若要求解該類問題，就必須用對應(yīng)于Stiff系統(tǒng)的指令

第三十五頁，共51頁。將值改成1000，ODE文件改寫成（vdp2.m）：

>>typevdp2.mfunctiondy=vdp2(t,y)mu=1000;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];第三十六頁，共51頁。選用專門用來解決Stiff系統(tǒng)問題的ODE指令，例如ode15s，來求解此系統(tǒng)并作圖顯示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)'); 第三十七頁，共51頁。由上圖可知，的變化相當(dāng)劇烈（超過），這也就是Stiff系統(tǒng)的特點(diǎn)第三十八頁，共51頁。繪制二維平面相位圖：ode15s02.m若要產(chǎn)生在某些特定時間點(diǎn)的狀態(tài)變量值，則調(diào)用ODE指令的格式可改成：

[t,y]=solver('odeFile',[t0,t1,…,tn],y0)其中[t0,t1,…,tn]即是特定時間點(diǎn)所形成的向量[t,y]=ode15s('vdp2',[03000],[21]');subplot(1,1,1);plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)')第三十九頁，共51頁。第四十頁，共51頁。ODE指令的選項(xiàng)ODE指令可以接受第四個輸入變量，代表積分過程用到的各種選項(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…

第四十一頁，共51頁。也可以只改變一個現(xiàn)有的options中某欄位的值，格式如下：newOptions=odeset(options,'name',value);若要獲取某欄位的值，可用odeget，格式如下：

value=odeget(otpions,'name');當(dāng)使用odeset指令時，若無輸入任何變量名，則odeset會顯示所有變量名及其值，并以大括號代表預(yù)設(shè)值第四十二頁，共51頁。>>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]第四十三頁，共51頁。由odeset產(chǎn)生的ODE選項(xiàng)

類別欄位名稱資料形態(tài)預(yù)設(shè)值說明誤差容忍度之相關(guān)欄位RelTol正數(shù)量相對誤差容忍度AbsTo1正數(shù)量或向量絕對誤差容忍度積分輸出之相關(guān)欄位OutPutFcn字串‘odeplot’輸出函式（若ODE指令無輸出變量，則在數(shù)值積分執(zhí)行完畢后，MATLAB會呼叫此輸出函數(shù)）OutputSel索引向量全部ODE指令之輸出變量的索引值，以決定那些輸出變量之元素將被送到輸出函式Refine正整數(shù)1或4(forode45)Refine=2可產(chǎn)生兩倍數(shù)量的輸出點(diǎn)，Refine=3可產(chǎn)生三倍數(shù)量的輸出點(diǎn)，依此類推Statson或offoffSt