常微分方程初值問題的MATLAB解法

1、用Matlab求常微分方程（ODE）的初值問題（IVP）本節(jié)考慮一階常微分方程u=f(t,u)u(t0)=U0t0<t<T(1.1)的數(shù)值求解問題，包括算法公式及編程問題。對(duì)一階常微分方程組的初值問題U'=fi(t,Ui,U2,，Um)U；=f2(t,Ui,U2,，Um)1Um二fm(t,Ui,U2,Um)Q。)=u°+/丁U2(t0)=U2t。二t_TJm(t0)=Um(1.2)可以通過引入列向量U,U0,f化成類似（1.1）的形式(1.3)U=f(t,U)t0<t_TiU(t0)=u0其中U(t)=“1(t)、U2(t)*,U0=%1收)、U2(t0)*

2、,f(t,U)=f1(t,U1,U2,，Um)"f2(t,U1,U2,，Um)I-(1.4)<Um(t)/um(t0)/fm(t,U1，U2,，Um)/另外一個(gè)高階常微分方程的初值問題(1.5)U(m)=f(t,U,U,U，U(m)t0：t<T(m)(m-1)U(t0)-U0,U(t0)-U0,U(t0)-U0也可以通過變換U1=U,U2=U：U3=U二，Um=U（m,）化成維微分方程組:U1二U24=U3«（1.6）Um4UmUm=f（t,U1，U2，,Um）我們?cè)谠O(shè)計(jì)算法時(shí)通常先對(duì)一維一階常微分方程（1.1）進(jìn)行，然后再將這個(gè)算法寫成適合求解（1.3）的向量

3、形式,并以向量形式來進(jìn)行編程。1、非剛性系統(tǒng)與剛性系統(tǒng)當(dāng)f（t,U）=AU+G（t）時(shí)微分方程組（1.3）變成U=Au（t）（1.7）如果系數(shù)矩陣AWRm2特征值九1,九21，九m滿足同之|%|2km|,對(duì)應(yīng)的特征向量為）,，Vm，則（1Q具有通解(1.8)m.4u(t)=xCjeVj(t)ji其中®(t)為(1.7)的一個(gè)特解。如果(1.7)中的矩陣滿足(1.9)(1) Re(j)<0j=1,2,m(2) s=Re(i)|.|Re(m).1就稱微分方程組(1.7)是剛性的或壞條件的或病態(tài)的°s稱為剛性比。對(duì)于一般的一階常微分方程組(1.3),如果多元向量值函數(shù)f(t

4、,u)又寸(1.4)中向量u的分量的Jacobi矩陣甘1甘1.生、cu1cu2dm況戔況J(t)=CUiCU2自m-s+勖S%編正<CUiCU2um/uRt)(1.10)的特征值X1(t),%(t),，九m(t)對(duì)于求解區(qū)間上的任何twt0,T滿足(1.9)式，則稱微分方程組(1.3)為剛性的。2、解法器及調(diào)用格式解法器適用類型何時(shí)使用使用算法ode45NonstiffMostofthetime.Thisshouldbethefirstsolveryoutry.Runge-Kutta(4,5)格式ode23NonstiffIfusingcrudeerrortolerancesorsolv

5、ingmoderatelystiffproblems.Runge-Kutta(2,3)格式ode113NonstiffIfusingstringenterrortolerancesorsolvingacomputationallyintensiveODEfile.Adams-Bashforth-MoultonPECE格式ode15sStiffIfode45isslowbecausetheproblemisstiff.ode23sStiffIfusingcrudeerrortolerancestosolvestiffsystemsandthemassmatrixisconstant.ode23t

6、ModeratelyStiffIftheproblemisonlymoderatelystiffandyouneedasolutionwithoutnumericaldamping.ode23tbStiffIfusingcrudeerrortolerancestosolvestiffsystems.有如下三種調(diào)用方法，其中soker是上面三個(gè)解法器的名子。(1) t,y=solver(odefun,tspan,y0);(2) t,y=solver(odefun,tspan,y0,options);(3) t,y=solver(odefun,tspan,y0,options,p1,p2.);它們

7、輸入?yún)?shù)是：(1) odefun:是一個(gè)函數(shù)句柄，指向初值問題(1.3)中的函數(shù)f(t,U),它的基本形式是dydt=f(t,y),其中dydt,y是列向量，而t是標(biāo)量。(2) tspan：是一個(gè)向量，用于指定求解區(qū)間tspan(1),tspan(end)。要求這個(gè)向量的分量是有序的，或者單調(diào)上升或者單調(diào)下降。y0:是初值問題(1.3)中的初始列向量同。(4)options：用于對(duì)解法器缺省的求解參數(shù)進(jìn)行設(shè)置。如果不想改變?nèi)笔≈?，可以用空矩陣?#39;來代替。MATLAB的ODE解法器設(shè)計(jì)了一個(gè)結(jié)構(gòu)變量用于設(shè)置解法器的各項(xiàng)公共參數(shù)，比如精度，步長等。其中已經(jīng)定義了一些缺省值，用無參數(shù)的ode

8、set命令可以列出所有的公共參數(shù)及其缺省值：AbsTol:positivescalarorvector1e-6絕對(duì)誤差RelTol:positivescalar1e-3相對(duì)誤差NormControl:on|off是否用向量范數(shù)控制誤差，缺省是按每個(gè)分量進(jìn)行誤差控制。OutputFcn:functionOutputSel:vectorofintegersRefine:positiveintegerStats:on|offInitialStep:positivescalar初始步長MaxStep:positivescalar最大步長，缺省為求解區(qū)間的十分之一。BDF:on|off是否利用向后差分來

9、求Jacobi矩陣。MaxOrder:1|2|3|4|5Jacobian:matrix|function對(duì)剛性方程組的解法器可以提供相應(yīng)的Jacobi矩陣信息JPattern:sparsematrixVectorized:on|offMass:matrix|functionMStateDependence:none|weak|strongMvPattern:sparsematrixMassSingular:yes|no|maybeInitialSlope:vectorEvents:function每個(gè)解法器可以設(shè)置的參數(shù)如下表：Parametersode45ode23ode113ode15so

10、de23sode23tode23tbRelTol,AbsTol,NormControlqq4qV4OutputFcn,OutputSel,Refine,StatsVV4EventsV_qV4MaxStep,InitialStepVqV4Jacobian,JPattern,Vectorized-7V4Mass77JMStateDependenceJJ-JMvPattern-MassSingular-InitialSlope-MaxOrder,BDF-(5)輸入?yún)?shù)p1,p2.等等是可選的參數(shù),并且由解法器彳遞給函數(shù)odefun.輸出參數(shù)的含義是：(1)t:節(jié)點(diǎn)列向量。(2)y：在節(jié)點(diǎn)處的解的近似

11、值數(shù)組，第i行對(duì)應(yīng)于t的第i行那個(gè)節(jié)點(diǎn)處的近似解。即輸出是如下形式:節(jié)點(diǎn)ty1(t)y2(t)ym(t)例題1、分別用步長h=0.05,0.1,0.125,0.2來求初值問題在區(qū)間0,8上的近似解，并與準(zhǔn)確解u=Ji+2t進(jìn)行圖形上的比較。2tu=u一下u(0)=10：t<8(1.11)解：利用如下代碼，注意其中作為輸入?yún)?shù)的函數(shù)的代碼的編寫方法，以及ODE解法器ode45的調(diào)用方法。functionzzzh=0.1;y0=1;tspan=0:h:8;t,y=ode45(f,tspan,y0);t,yholdon;plot(tspan,y,'-k',tspan,fu(ts

12、pan),'-b');title('SolutionofIVPineg1');xlabel('timet');ylabel('solutionu');legend('numericalsolutons','truesolutions');functiondudt=f(t,u)dudt=u-2*t./u;functionu=fu(t)u=sqrt(1+2*t);由于結(jié)果較多我們只畫出它們的圖形如圖1。從圖形中可以明顯看出，當(dāng)時(shí)間變量較大時(shí)比如tA5所求的近似解嚴(yán)重偏離了真解。這說明初值問題(1.11)

13、在時(shí)刻0有較好的穩(wěn)定性，但是在時(shí)刻4以后可能穩(wěn)定性就不太好了，為些我們利用(1.11)構(gòu)造一個(gè)初值問題：(1.12)u'=u-胃4<tE12u(4)=3它與初值問題(1.11)有相同的解。將上述代碼適當(dāng)改動(dòng)來求解初值問題(1.12),并畫出圖形如圖2,發(fā)現(xiàn)還是在初始時(shí)刻4附近有較好的穩(wěn)定性，所以我們?cè)谇蠼鈺r(shí)不能希望在很大的區(qū)間上求解。1614Sow2nmlvPii>egl一numericaladhilonLnjesoluiionfi圖1圖2例2、求解二階VanderPol微分方程所構(gòu)成的初值問題在區(qū)間0,20上的近似解。其中參數(shù)N是一個(gè)正數(shù),當(dāng)較大時(shí)這個(gè)方程是剛性的。(1.

14、13)uJ.；(1u2)uu=00<t<20u(0)=2,u(0)=0解：這是一個(gè)高階方程，要用變量代換化成向量形式的一階方程組形式(1.3),令v=u'則化成方程組rfu=v2vN(1u)u+u=00<t<20(1.14)u(0)=2,v(0)=0再令u(t)=L(t),v(t)T,f(t,1)=v,N(1u2)vuT則(1.14)的一階向量形式為(1.3)。下列代碼中編寫了向量值函數(shù)f(t,u),由還有一個(gè)可選參數(shù)N故要注意解法器的調(diào)用形式。輸出矩陣u的每一行代表解向量在每個(gè)節(jié)點(diǎn)處的近似函數(shù)值。對(duì)于本題來說輸出矩陣u的第二列相當(dāng)于方程(1.13)解函數(shù)u(t

15、)導(dǎo)數(shù)的近似值?？梢詫?duì)N=0.5,1,2,5,10,100,1000進(jìn)行試驗(yàn)解法器ode45的求解結(jié)果及圖像。functionzzz%ode13.mmu=10;mustr=num2str(mu);h=0.1;u0=2,0'tspan=0:h:20;t,u=ode45(f,tspan,u0,mu);t,ufigure;holdon;plot(t,u(:,1),'-k',t,u(:,2),'-b');title(strcat('例2中的初值問題mu=',mustr);xlabel('timet');ylabel('so

16、lutionu');legend('數(shù)值解,'導(dǎo)函數(shù)數(shù)值解');holdoff;functiondudt=f(t,u,mu)dudt=u(2),mu*(1-u(1)八2)*u(2)-u(1)'倒3中的初值問題產(chǎn)0.1246B1012141E1B20time125150552-5tO.。-1.“'2.=8虧七|3教值解身函救第值解口-1-2-3n8虧七6后臺(tái)PHJ-32表1:當(dāng)mu=0.5時(shí)的幾個(gè)解節(jié)點(diǎn)tu(t)u'(t)02.000000.10001.9905-0.18550.20001.9638-0.34450.30001.9223-0

17、.48140.40001.8681-0.60080.50001.8026-0.70680.60001.7271-0.80340.70001.6422-0.89400.80001.5484-0.98130.90001.4459-1.06771.00001.3349-1.15461.10001.2150-1.24431.20001.0858-1.3386如果再增大N,則ode45解法將無法求出較好的數(shù)值解，原因是當(dāng).斯11y11102661012141610301i1N較大時(shí)方程組(1.14)的Jacobi矩陣模最大特征值與模最小特征比值比較大，或者說剛性比大。(1.15)f_f01、£

18、u一2Nuv-1N(1-u2)，牙001在時(shí)刻0處，當(dāng)N=100時(shí)(1.15)為：一=,其特征值為：-0.0033,-299.9967,剛性比約為cut=01-300,90909,非常大，這種方程不能用顯式的四階RK方法來求解，必須用某種隱式方法才能求解。下面代碼是利用ode23s解法器來求解VanderPol方程當(dāng)N=10,100,1000,10000時(shí)的解，由于剛性比較大，例工中的初值同超K=10W對(duì)應(yīng)絕對(duì)值較小特征值的解分是變化較緩慢，故應(yīng)當(dāng)采用很大的求解區(qū)間，才能觀察到解的變化情況，因此取求解區(qū)間為0,5用。同學(xué)們自己也可以使用較小的區(qū)間來試驗(yàn)，看看數(shù)值解的模擬圖像能否得到真解的性態(tài)。

19、由于數(shù)據(jù)比較多，故只給出數(shù)值解模擬圖像。functionzzzmu=100;mustr=num2str(mu);h=0.1;u0=2,0'tspan=0:h:5*mu;t,u=ode15s(f,tspan,u0,mu);t,u;figure;holdon;plot(t,u(:,1),'-b');title(strcat('例2中的初值問題mu=',mustr);xlabel('timet');ylabel('solutionu');legend('數(shù)值解')；holdoff;functiondudt=f(t

20、,u,mu)dudt=u(2),mu*(1-u(1)A2)*u(2)-u(1)'例3、求解如下非剛性初值問題，它可以描述剛性物體不受外力作用時(shí)的運(yùn)動(dòng)狀態(tài)。y1=y2y3y1(0)=0y2=-y1y30；t_12y2(0)=1(1.16)3=0.51丫2乂)3(0)=1解：由于是非剛性問題故可以使用ode45來求解，這個(gè)例子說明怎樣設(shè)置解法器的options選項(xiàng)中的參數(shù)。具體代碼如下：functionzzzoptions=odeset('RelTol',1e-4,'AbsTol',1e-41e-41e-5);tspan=0,12;y0=0,1,1;T,Y=

21、ode45(rigid,012,011,options);T,Yplot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')title('例題3中的圖)；xlabel('時(shí)間t');ylabel('解y(t)');legend('數(shù)值解y_1(t)','數(shù)值解y_2(t)','數(shù)值解y_3(t)');functiondy=rigid(t,y)dy=zeros(3,1);%acolumnvectordy(1)=y(2)*y(3);d

22、y(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);解的數(shù)值模擬圖像為：練習(xí)1、求解例1中的vanderPol微分方程，其中N=0,1/16,1/4,1/2,1,2,4,10,20,100,1000,并到互聯(lián)網(wǎng)上去搜索這個(gè)微分方程是可以描述哪些問題?練習(xí)2、用上面的所有方法求解如下初值問題的解2.y-0.1(1-y)yy=0y(0)=1,y(0)=0x(t)=-=r,ry=-2r(3)y-2y3=0y(1)-y(1)-1,1：x<2y=2z3x(4)z=2yz.y(0)=z(0)=1”,x(0)=0.4,x(0)=0y(0)=0,y(0)=2練習(xí)3、求解例1中的vanderPol微分方程，并到互聯(lián)網(wǎng)上去搜索這個(gè)微分方程是可以描述哪些問題？其中:二二0,1/16,1/4,1/2,1,2,4,10,20,100,1000練習(xí)4、Lorenz吸引子是在渾沌世界中經(jīng)常用到的常微分方程組，它實(shí)際上描述了容器內(nèi)的液體被加熱后在容器內(nèi)的流動(dòng)情況，可以用于描述大氣的對(duì)流情況。x=p(y-x)y=rx-y-xz(1.17)z二xy