常微分方程的數(shù)值解法分析解析

1、第八章常微分方程的數(shù)值解法i內(nèi)容要點考慮一階常微分方程初值問題：dxf (x, y)y(xo)=yo微分方程的數(shù)值解：設微分方程的解y(x)的存在區(qū)間是a,b,在a,b內(nèi)取一系列節(jié)點a=xo<xi<<xn=其中hk=xk+i-xk；(一般采用等距節(jié)點，h=(b-a)/n稱為步長)。在每個節(jié)點xk求解函數(shù)y(x)的近似值：y-y(xk),這樣yo,yi,.,yn稱為微分方程的數(shù)值解。用數(shù)值方法，求得f(xk)的近似值y.再用插值或擬合方法就求得y(x)的近似函數(shù)。(一)常微分方程處置問題解得存在唯一性定理對于常微分方程初值問題:dXH(x,y)如果：y(xo).。(1)在x。，

2、人,卜%的矩形內(nèi)f(x,y)是-個二元連續(xù)函數(shù)。(2)f(x,y)對于y滿足利普希茨條件，即f(x,yi)»(x,y2)|L|yi21則在x。|x|c上方程.£=f(x,y)的解存在且唯-,這y(xo)3y0里C=min(A-xo),xo+B/L),L是利普希茨常數(shù)。定義：任何一個一步方法可以寫為y|Jyk(xk,yk,h),其中同(xk,yk,h)稱為算法的增量函數(shù)。收斂性定理：若一步方法滿足：是p解的.(2)增量函數(shù)圖x,y.h)對于y滿足利普希茨條件.(3)初始值yo是精確的。則|y(kh)，y(x)倡O(hp),kh=x-xo，也就是有hlim，k.y(x)號kh-

3、x-xQ(一)、主要算法1.局部截斷誤差局部截斷誤差：當y(xk)是精確解時，由y(xk)按照數(shù)值方法計算出來的丫息的誤差y(xk+i)-yk稱為局部截斷誤差。注意：yk+i和y屆i區(qū)別。因而局部截斷誤差與誤差ek+i=y(xk+i)-yk+i不同。如果局部截斷誤差是O(hp+i),我們就說該數(shù)值方法具有p階精度1 .顯式歐拉格式:，k 1=yk hf (xk, yjk=1,2, n-1.y(xo) = yo顯式歐拉格式的特點：(1)、單步方法；(2)、顯式格式；(3)、局部截斷誤差O(h2)因而是一階精度隱式歐拉格式y(tǒng)k hf(Xk i,yk i) y(%) = y。2 .兩步歐拉格式:&#

4、39;k 1 =%2hf(Xk，yk)k=1,2, n-1.y(X。) = y。兩步歐拉格式的局部截斷誤差 o(h3),因而是二階精度3 .梯形方法：yk 1,ykyk h f (Xk 13.改進的歐拉方法:預測值：yky(xo) = yoykh f (Xk, yk)校正值:ykykhyp BykhDf(Xk,yk)f(Xk i,yk i) . f (Xk,yk)或改寫為ICBykhBf(xlllyp)1.、yki二2(ypyJ4、梯形方法與改進歐拉方法的截斷誤差是O(h3),具有二階精度。5、龍格-庫塔法的思想1).二階龍格-庫塔法計算公式:K1. (Xk, yk)心(XkHph, yk 3

5、phK1)32yki時，得一簇龍格-庫塔公式，其局部截斷誤差均為O(h3)都是二階精度。特別2取廈1,P.1，就是改進歐拉公式。2取p=1,得二階龍格一庫塔法為:Ki=f(x.yk)K2Bf(xkBh,ykHhK1),稱為二階中點格式。ykhK22)、經(jīng)典龍格-庫塔格式(也稱為標準龍格-庫塔方法)：Ki=f(Xk,yk)hhK2=f(Xk,ykKi)22K3=f(Xkh,ykh.)22K4=f(xh,yh.)ykh(K12K22K3K4)6四階龍格一庫塔方法的截斷誤差為h5I,具有四階精度。值問題用四階龍格-庫塔方法計算，其精度均滿足實際問題精度要求般一階常微分方程初3) .變步長龍格庫塔方法

6、：從節(jié)點Xk出發(fā)，以步長h據(jù)四階龍格庫塔方法求出一個近似值ykj|然后以步長h/2求出一個近似值y|,得誤差事后估計式：(h/2)1(h/2)(h)yk-yk1'.(yki-丫)根yllliy|來選取步長h。4) .RKF格式：變步長龍格庫塔方法，因頻繁加倍或折半步長會浪費計算量。Felhberg改進了傳統(tǒng)龍格-庫塔方法，得到RKF格式，較好解決了步長的確定，而且提高了精度與穩(wěn)定性，為Matlab等許多數(shù)值計算軟件采用。4/5階RKFB式：由4階龍格-庫塔方法與5階龍格-庫塔方法結合而成。25216K11408K321972565410116135ayy6656.1282528561K

7、3K56430K4-555J50K1-f(xn,yn)K2=f(XnLyn"1)44K3=f(xn8h，yn329K132K2)12h1932K4=f(XnITyn-7200K12197K2K2197K5=f(Xnh,yn4393680K18K2巫0216513f(xnh,yn-_8Kl2K2-3544K322725658454104K3)K61859410411K4-40K4)1/4為精度要求,y二階顯式Adams方法:三階顯式Adams方法:四階顯式Adams方法:hyk 1 = yk 23fk - fk,h .yk 1 = yk 12 23 fk-16fk5fk/;yk 1 =

8、 y2h4-9 fkj3 37 fk/ 一59 fk55 fk.Felhberg得到的最佳步長hs,其中h為當前步長,若s<0.75,步長折半；若s>1.5步長加倍6.亞當斯方法(Adams)1) .顯式AdamspJ法:記：fkHf(xk,yk);三階隱式Adams方法:yk 23fk - fk 1;yk .5fk 1 8fk - fk四階隱式Adams方法:yk 1 = y24(fk -5f19fk 9fk1)2) .隱式Adams方法二階隱式Adams方法:3) .Adams預報-校正系統(tǒng)：先用顯式格式作為預測值，再用隱式格式來校正。預測化yk1=yk，(55fk-59fk3

9、7f_-9f-)校正值:ykLyk24(f-5f19fk94).改進的Adams預報-校正系統(tǒng):改進：m預測：Pk校正值:改進：yk1(ckii)0,(同時fi一Pki)時趨向于準確解y(x.，(55fk.59fk437fy_9f.),m251,、(ck-pk)270Ck1=y,9fxk119ck12707.收斂與穩(wěn)定性對于固定的xk.x0Bkh,如果數(shù)值解yk當19fk-5f則稱該數(shù)值方法方法是收斂的。ym(m>k)上產(chǎn)生的如果一種數(shù)值方法，在節(jié)點值yk上大小為事擾動，于以后各節(jié)點值偏差均不超過、則稱數(shù)值方法是穩(wěn)定的.8.剛性方程組：考慮n階常微分方程組:)n的實部Re(i)<0

10、,i=1,2,n.,ode23、ode11& ode15般，隱式格式較顯式格式有較好的穩(wěn)定性。若矩陣A的特征值的剛性比，當剛性比很大時稱*式為剛性方程組。剛性方程組需采用穩(wěn)定性較好的方法。二.MATLAB的有關命令t,x=solver('f',ts,x0,optiOnst：輸出的自變量值,x:輸出的函數(shù)值；solver:包才5ode4ssode23s.非剛性系統(tǒng):ode23:用2階(及3階)龍格-庫塔算法。ode45:用4階(及5階)龍格-庫塔-算法。ode113(Adams-Bashforth-MoultonPECE侈步方法。剛性系統(tǒng)：ode15s(GearTj法)多

11、步方法ode23sL階modifiedRosenbroackformula并步。ode23t(trapezoidalrule)solveDAEs.ode23tb(TR-BDF2)低精度算法。f：由待解方程寫成的m-文件名ts=t0,tf,t0、tf為自變量的初值和終值x0:函數(shù)的初值options=odeset(reltol',rt,')ab1sltoat：'角制為設定的相對誤差和絕對誤差.(缺省時設定相對誤差10-3,絕對誤差10-6)一、/A-汪忠：1、在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-文件中的待解方程組應以x的分量形式寫成.2、使用Matlab

12、軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.四.重點算法Matlab程序Euler格式、functionx,y=naeuler(dyfun,xspan,y0,h)x=xspan:h:xspan(2);y(1)=y0;forn=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);endx=x'y=y'隱式Euler格式functionx,y=naeulerb(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1y(n+1)=iter(dyfu

13、n,x(n+1),y(n),h);endx=x'y=y'functiony=iter(dyfun,x,y,h)y0=y;e=1e-4;K=1e+4;y=y+h*feval(dyfun,x,y);y1=y+2*e;k=1;whileabs(y-y1)>eyi=y;y=y0+h*feval(dyfun,x,y);k=k+1;ifk>K,error('迭代發(fā)散');endend改進Euler格式functionx,y=naeulerg(dyfun,xspan,y0,h)x=xspan:h:xspan；y(1)=y0;forn=1:length(x)-1K1

14、=feval(dyfun,x(n),y(n);y(n+1)=y(n)+h*K1;K2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(K1+K2)/2;endx=x'y=y'4階經(jīng)典Runge-Kutta格式functionx,y=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1K1=feval(dyfun,x(n),y(n);K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1);K3=feval(dyfun,x(n)+h/2,y(n

15、)+h/2*K2);K4=feval(dyfun,x(n+1),y(n)+h*K3);y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6;endx=x'y=y'變步長4階經(jīng)典Runge-Kutta格式functionx,y=nark4V(dyfun,xspan,y0,e,h)ifnargin<5,h=(xspan(2)-xspan(1)/10;endn=1;x(n)=xspan(1);y(n)=y0;y1,y2=comput(dyfun,x(n),y(n),h);whilex(n)<xspan(2)-eps;ifabs(y2-y1)/10>ew

16、hileabs(y2-y1)/10>eh=h/2;y1,y2=comput(dyfun,x(n),y(n),h);endelsewhileabs(y2-y1)/10<=eh=2*h;y1,y2=comput(dyfun,x(n),y(n),h);endh=h/2;h=min(h,xspan(2)-x(n);y1,y2=comput(dyfun,x(n),y(n),h);endn=n+1;x(n)=x(n-1)+h;y(n)=y2;y1,y2=comput(dyfun,x(n),y(n),h);endx=x'y=y'functiony1,y2=comput(dyfun

17、,x,y,h);y1=rk4(dyfun,x,y,h);y21=rk4(dyfun,x,y,h/2);y2=rk4(dyfun,x+h/2,y21,h/2);functiony=rk4(dyfun,x,y,h);K1=feval(dyfun,x,y);K2=feval(dyfun,x+h/2,y+h/2*K1);K3=feval(dyfun,x+h/2,y+h/2*K2);K4=feval(dyfun,x+h,y+h*K3);y=y+h*(K1+2*K2+2*K3+K4)/6;五.試驗例題例1不同方法的精度比較用多種方法解下述初值問題，并與其準確解ylIBx進行比較。y=-yx1,0：x<

18、;0.6y（0）=1解：取步長h=0.1,xk=kh（k=0,1，,6）。各方法計算結果及絕對誤差見表（1）、（2）、（3）表（1）xn歐拉公式改進的歐拉公式格階標準龍格一庫塔公式y(tǒng)n誤差yn誤差yn誤差。1.0000000101011.00000043KL尸1.00500001.6xlCH1.004837500121.010000&7x1戶1.0190252.9x101.018730901.5x100131.0290001.2kIO-31.0412184.0x101.0408184220IX10-741.0561001.4xlO-31.07080243x11.070320292.4x

19、1a-70151.0904901.6xlO-51.1070765.5x101.106530932.7x10-7IH1.81.1494045-1.14881193，力xltH表（2）四階阿當姆斯公式nXn顯示公式隱式公式y(tǒng)n誤差yn誤差130.3取自準確解一1.040818011x1尸40.41.070322921.070319663.9x10-?50.51.106535481.106530415.2x1060.61.148814811.1488110165x10備注y1,y2,y3取自準確解y1,y2取本題/,”I關于y為線性，代入隱式公式，可解出yn+i,因此可直接求隱式公式之解。表(3)四

20、階阿當姆斯預測一校正公式nXn顯示公式隱式公式Yn誤差yn誤差40.41.070319921.3xl0-71.070320140.PX1O-750.51.106530273.510-71.106530771.1x1Q-760.61.148811036.1xIQ-71.148811751.2x10-7備注y1,y2,y3由四階標準龍格一庫塔公式提供對比以上各表數(shù)據(jù)知，在相同步長下求解同一問題時，方法階數(shù)越高，解的精度也越高，一階的歐拉公式精度最低，格階標準龍格一庫塔公式的精度大大高于改進的歐拉公式。四階阿當姆斯方法、顯式的精度略低于隱式。同是阿當姆斯預測一校正系統(tǒng)，帶誤差修正的精度又高于不帶誤差

21、修正的。從計算量看，四階龍格一庫塔法計算量最大，每前進一步要計算4次函數(shù)值f,而誤誤差修正的阿當姆斯預測一校正法與它精度差不多的，每前進一步只計算兩次函數(shù)值f,所以后者是可取的。但它是四步法，必須用其它方法提供出始值，程序略復雜些。例2(步長的計算結果的影響)用歐拉公式求下述初值問題在x=1處的近似解，并與準確解y(1)=1比較。J*-1000(y-#)+21yQ)=1解分別取步長h=10-1、h=1。-2、h=1,3、h=10-4計算，結果見下表。由表中數(shù)據(jù)可見，當卜=,時結果完全失真，而取,二10比力=10計算量增加了十倍，但解的精度卻基本一樣，可見取。二1CH太浪費計算量了。Xhy(1)

22、的計算值io-10904382075CBxWw1尸>10翔0.9599993000011戶0999999?00000上述結果差異很大的原因在于歐拉公式的絕對穩(wěn)定區(qū)間為（-2, 0）步長h應滿足眉(。),對本題，f (x,y)B 0 0OyBx2)g2x,故應取h滿足-2 <-1000 <0即0.，<2x10，?？梢娙×?=1。時歐拉公式是數(shù)值不穩(wěn)定的，導致結果失真，而取和= 1。,都滿足穩(wěn)定性要求，可用于求解。此例說明，求解微分方程數(shù)值解一定選取注意步長，過大則導致解的失真，過小又 會使計算量大增。究竟取多大步長才合適，不僅取決于所采用的數(shù)值方法，還決定于待 解微分方程

23、本身的特性。例3: （Matlab不同命令之差異）取h=0.2分別用不同的 Matlab命令解方程y - 2t/yy(1) = 1(0<t<4).其準確解為：y = % 1 2t解：>> t,y=ode45(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans =45.00000000000000 0.00020257808813>> t,y=ode23(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans

24、=13.00000000000000 0.19048360113013>> t,y=ode113(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans =18.00000000000000 0.00973277413962>> t,y=ode23t(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans =18.00000000000000 0.03919668151270>> t,y=ode23s(odefu

25、n,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans =81.00000000000000 2.54374270001824>> t,y=ode23tb(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans =15.00000000000000 0.24308870757312>> t,y=ode15s(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).A2/n);n,eans

26、 =22.00000000000000 0.45509043060378可見ode45精度較高，計算量較大，ode23計算量小，但誤差大，ode113適中，用剛性方程組揭發(fā)解非剛性問題不適合，特別ode23s計算量大且但誤差大。例4:解剛性方程組:-0.01y-99.99z,y(1)= -100z,z(0)=1其準確解為：y!JIII1xE19l0x,zjj0'就h=0.01、h=0.02使用改進歐拉格式及方程組，并在-張圖上繪出其圖形。functionx,y=naeuler2s(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length

27、(y0),length(x);y(:,1)=y0(:);forn=1:length(x)-1K1=feval(dyfun,x(n),y(:,n);y(:,n+1)=y(:,n)+h*K1;K2=feval(dyfun,x(n+1),y(:,n+1);y(:,n+1)=y(:,n)+h*(K1+K2)/2;endx=x'y=y'%用歐拉格式解常微分方程組的m文件functionf=dyfun(x,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:);%函數(shù)m文件clear;x,y=naeuler2s(dyfun,0500,2,1,0

28、.01);plot(x,y),axis(-50500-0.52)gtext('y(x)h=0.01'),gtext('z(x)h=0.01')%步長為0.01時解方程并畫圖clear;x,y=naeuler2s(dyfun,0500,2,1,0.02);holdonplot(x,y),gtext('y(x)h=0.02'),gtext('z(x)h=0.02')%步長為0.02時解方程并畫圖。clear;x=0:0.01:500;y=exp(-0.01*x)+exp(-100*x);z=exp(-100*x);holdonplot

29、(x,y,'r',x,z,'r'),gtext('y(x)真根'),gtext('z(x)真根')畫出真根圖形運行結果真根與h=0,01時結果重合，h=0.02時結果是錯誤的。故剛性方程組應采用穩(wěn)定性較好方法。例6:數(shù)值實驗(摘自蔡)觀察歐拉顯示方法的收斂性。內(nèi)容：取h=0.1,0.01.用歐拉顯式方法求解：用歐拉法解方程"|，3(0<x<1).y(1)=1計算到y(tǒng)(2),并與精確解丫(*)/酷2。2)/3_3,2比較。MATLAB程序t=1;y=1;h=0.1fork=1:10;y=y+h*t*y(1/3)

30、t=t+hend運行結果0.10.010.0010.0001y=(t2+2)/3)3/2(準確解)1111111.11.11.11751.10681.10681.106816606308381.21.21361.23931.22771.22791.227879593566201.31.34161.37621.36391.36411.364135990288361.41.48491.52941.51621.51651.516564538686041.51.64471.69981.68571.68611.686170601183731.61.82171.88831.87341.87391.8739

31、81856902571.72.0172.09622.08042.08102.081044689573001.82.23192.32432.30762.30832.308421169608421.92.46712.57392.55632.55712.557186539930162.02.72392.81772.82742.82832.82842712474619歐拉顯示方法是收斂的。例7：實驗目的：比較不同算法用于“非光滑”解時的結果內(nèi)容：用歐拉顯式方法和2階龍格，庫塔算法求解:|sinkx|, y(0)副，計算 y(福)并在一張圖上畫出解曲線k=6時的程序如下：cleart=0;y=1;h=p

32、i/60;fork=0.1:0.1:2;y=y+h*abs(sin(6*t)t=t+hplot(t,y,'o'),holdon,endt0=0;tf=1/3*pi;x0=1;k=pi/60t,x=ode23('ex545f,t0:k:tf,1)plot(t,x),holdon,fori=0:pi/60:pi/3ifi<=pi/6;y1=7/6-1/6*cos(6*i)elsey1=9/6+1/6*cos(6*i)endplot(i,y1,'*'),holdon,end例8實驗目的：觀察歐拉顯式方法的數(shù)值不穩(wěn)定性。內(nèi)容：取卜=0.05,0.04,0.

33、03,用歐拉顯式方法求解：B0y,y(0)Bo0,并畫出曲線，同樣用歐拉隱式方法重復求解上述問題。解：歐拉顯式方法程序略。因y10y是y的線性函數(shù)故用歐拉隱式方法可直接解出y。程序如下clear;h=0.01;t=0;y=100;fork=0:h:1;y=y/(1+50*h);t=t+h;plot(t,y,'b*'),holdon,endgtext('隱t=0.01')運行結果：100302010-1020-30 40-9046403530252015-rM1D353010100.20.40.E11.2-TO50a3DX)1004-=*4-*隱1=0.Dfl+:

34、+4-+4-+4-4-+4-LJ0.2D4D.G0.8121.2515020.JD.60.011.21J=-+-*-F號EJ03k如3Sno252D25t20»熄1<I.CG?*kumin”Q02O<06圖中標示為“隱00x”為隱式歐拉格式結果，其余是歐拉格式的運行結果，可見隱式歐拉格式較歐拉格式穩(wěn)定。2t2)(1 - cont) - -210-3，如果用長格式輸出y、本題解析解：y=(1MATLAB程序functionxdot=ex545f(t,x)u-2)；xdot=01;-10*x+01*u;%建立ex545f.m函數(shù)程序。clf,t0=0;tf=3*pi;x0t=

35、0;0;%給出初始值。t,x=ode23('ex545f,t0,tf,x0t)%調(diào)用MATLAB中的數(shù)值積分2階3階龍一庫公式。y=x(:,1);%丫是*的第2歹。forI=1:length(t);y2(I)=(1+2/(piA2)*(1-cos(t(I)-t(I)A2/(piA2)end%在數(shù)值積分的序列點上計算解析解的值。u=1-(t.A2)/(piA2);clf,plot(t,y,'-',t,u,'+',t,y2,'o')%繪出數(shù)值解、解析解的圖形。legend（'數(shù)值解輸入量解析解'）運行結果解析解與數(shù)值解看不出差

36、別，是因為本數(shù)值積分是按默認精度y2可看出它們的微小誤差。例10解微分方程組.y/gy2y3丫20y1y3ysmy1y2y1(0)=0羋(0)=1孚(0)=1解1、建立m-文件rigid.m如下：functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3)；dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0,tf=12,輸入命令：>>T,Y=ode45('rigid',012,011);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y

37、(:,3),'+')3、結果如圖.【例11】采用最簡潔格式的ODE文件和解算指令，研究圍繞地球旋轉的衛(wèi)星軌道。(1)問題的形成軌道上運動的衛(wèi)星，在Newton第二定律2drF=ma=m2-,和萬有引力定律dt2mM E3±r作用下，有 rd2r Mxyax-GMe3ay-GMe3,而rrx2臥2這里G等.672010J|(NM2/kg2)是引力常數(shù)，MEF.97即024(kg)是地球(2)構成一階微分方程組y24-GMeyi初始條件為Y(0)(0)y3y4yy3y4yi2 22.3/2(x y )y2(x2臥2嚴司0 0 Vy(0)E(5.14.2.1-5)(5.14

38、.2.1-6)的質(zhì)量。又假定衛(wèi)星以初速度Vy(0)H4000(m/s)在x(0)J通07(m)處進入軌道。平旱軌道圖(1)根據(jù)式(5.14.2.1-5)編寫最簡潔格式的ODE文件DYDt2.mfunctionYd=DYDt2(t,Y,flag,G,ME)%flag按ODEC件格式規(guī)定，必須是第三輸入宗量。對它的賦值由ode45指令自動產(chǎn)生。第4、5宗量是被傳遞的參數(shù)switchflagcase''%按規(guī)定：這里必須是空串。在此為"真"時，完成以下導數(shù)計算。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.A2);Yd=V;-G*ME*X/rA3;

39、otherwiseerror('Unknownflag'''flag'''.');end(2)對微分方程進行解算G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;%指定解算微分方程時間區(qū)間Y0=x0;0;0;vy0;%按式(5.14.2.1-6)給定初值向量也YY=ode45('DYDt2',tspan,Y0,口,G,ME);%第4宗量取缺省設置，第5、6宗量是被傳遞參數(shù)X=YY(:,1);%輸出Y第一列是位移數(shù)據(jù)x(t

40、)Y=YY(:,2);%輸出Y第二列是位移數(shù)據(jù)y(t)plot(X,Y,'b','Linewidth',2);holdonaxis('image')%保證x、y軸等長刻度，且坐標框恰包容圖形XE,YE,ZE=sphere(10);%產(chǎn)生單位球面數(shù)據(jù)RE=0.64e7;%地球半徑XE=RE*XE;YE=RE*YE;ZE=0*ZE;%坐標紙上的地球平面數(shù)據(jù)mesh(XE,YE,ZE),holdoff%繪地球示意圖【*例12帶事件設置的ODE文件及主程序編寫演示。本例將以較高精度計算衛(wèi)星經(jīng)過近地點和遠地點的時間，并在圖上標志。(1)ODE文件的編寫DY

41、Dt3.mfunctionvarargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)%DYDt3.m供主程序調(diào)用的ODE!數(shù)文件%本文件自帶三個子函數(shù)：f,fi,fev。%t,Y分別是自變量和一階函數(shù)向量，是最基本的輸入宗量。%flag第三輸入宗量，它專供解算指令(如ode45)作調(diào)用通知。%在運行中，解算指令會根據(jù)需要向flag發(fā)不同的字符串。%varargout是"變"輸出宗量。它由變維的元胞數(shù)組構成。每個元胞中可以存放指令%所產(chǎn)生的任意形式的數(shù)據(jù)。switchflagcase''%必須用空串符。情況為"真"時，計

42、算導數(shù)dY/dt=f(t,Y)。varargout1=f(t,Y,G,ME);%輸出為一個元胞，容納f子函數(shù)的一個輸出Ydcase'init'%必須用init'，情況為"真"時，傳送計算區(qū)間、初值、設置參數(shù)。varargout1:3=fi(tspan,Y0);%諭出為三個元胞，容納fi子函數(shù)的三個輸出量case'events'%必須用events'，情況為"真"時，設置事件性質(zhì)。varargout1:3=fev(t,Y,Y0);%諭出三個元胞，容納fev子函數(shù)的三個輸出量otherwiseerror(

43、9;Unknownflag'''flag'''.');end%functionYd=f(t,Y,G,ME)%十算導數(shù)子函數(shù)，被"父"函數(shù)DYDt耐用。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.A2);Yd=V;-G*ME*X/rA3;%functionts,y0,options=fi(tspan,Y0)煙置時間區(qū)間、初值、算法參數(shù)子函數(shù)，被"父"函數(shù)DYDt3調(diào)用。ts=tspan;y0=Y0;options=odeset('Events','on'

44、;,'Reltol',1e-5,'Abstol',1e-4);%F動。de45的"事件判斷"功能，設置相對誤差和絕對誤差。%functionvalue,isterminal,direction=fev(t,Y,Y0)%M牛判斷子函數(shù)，被"父"函數(shù)DYDt耐用。dDSQdt=2*(Y(1:2)-Y0(1:2)'*Y(3:4);%dDSQ慳離初始點的二乘距離u的時間導數(shù)du/dt,而u=(x-x0)A2+(y-y0)A2。value=dDSQdt;dDSQdt;%定義兩個穿越0的事件direction=1;-1;蹴一事

45、件：以漸增方式穿越0。第二事件：以漸減方式穿越0isterminal=1;0;蹴一事件發(fā)生后，終止計算；而第二事件發(fā)生后，繼續(xù)計算。(2)運行以下主程序G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;Y0=x0;0;0;vy0;t,YY,Te,Ye,Ie=ode45('DYDt3',G,ME,tspan,Y0);%<3>X=YY(:,1);Y=YY(:,2);plot(X,Y,'b','Linewidth',2);holdontext(0

46、,6e7,'軌道','Color','b')%產(chǎn)生藍色文字注釋axis('image');%保證x、y軸等長刻度，且坐標框恰包容圖形在三個事件發(fā)生點上畫標記plot(Ye(1,1),0.4e7+Ye(1,2),'rA','MarkerSize',10)plot(Ye(2,1),0.4e7+Ye(2,2),'bv','MarkerSize',10)plot(Ye(3,1),-0.4e7+Ye(3,2),'bA','MarkerSize',

47、10)咐巴軌道的半周期和全周期標在圖上text(0.8*Ye(3,1),-2e7+Ye(3,2),'t3='sprintf('%6.0f',Te(3)text(0.8*Ye(2,1),1.5e7+Ye(2,2),'t2='sprintf('%6.0f',Te(2)%在x-y坐標上畫地球XE,YE,ZE=sphere(10);RE=0.64e7;XE=RE*XE;丫E=RE*YE;ZE=0*ZE;mesh(XE,YE,ZE)text(1e7,1e7,'地球','Color','r'),

48、holdoff%產(chǎn)生紅色文字注釋例13范例：狀態(tài)轉移方程組模型計算機網(wǎng)絡可靠性分析問題隨著計算機通信網(wǎng)絡系統(tǒng)特別是Internet網(wǎng)絡日益廣泛的應用，提高系統(tǒng)的可靠性意義重大.研究和分析具有實用性的高可靠計算機通信網(wǎng)絡系統(tǒng)，是國際上非?；钴S的一個研究方向.(1)問題及假設我們將研究的計算機通信網(wǎng)絡系統(tǒng)稱為無冗余的防火墻協(xié)議系統(tǒng).計算機隨時可能發(fā)生三個事件一無故障、間歇故障和永久故障.因此，計算機一般處于三種工作狀態(tài)：無故障工作、帶故障工作和不工作.這三種狀態(tài)之間的轉移過程如圖3.9所示.要求建立該系統(tǒng)的狀態(tài)轉移模型，并進行可靠性分析.(2)分析與模型圖3, 9無冗余的防火墻協(xié)議系統(tǒng)狀態(tài)轉移討程

49、該問題屬于狀態(tài)轉移問題，利用馬爾科夫狀態(tài)轉移原理，用RO),P2(t)和)分別表示系統(tǒng)處于無故障工作、帶故障工作和不工作三種狀態(tài)的概率，則有以下狀態(tài)轉移方程組：(即2)p(t)碣同園)P2(t)'2(t)p t/2)Pi(t)初始條件P(0),P2(0),P3(0).0,°,參數(shù)取值. ：10" 10"|：104 10： -O.1這是一個帶參數(shù)的微分方程組模型.dP2(t)(t)，2(t)(3)模型求解求解這個帶參數(shù)的微分方程組模型有兩種方法，一是用特征根法求解析解，另一個是用數(shù)值解法求數(shù)值解.分別求解如下：假定模型中的參數(shù)取下限值，即"=10；

50、10一4,=0.01.解析解法利用常系數(shù)線性微分方程組的特征根法.MATLAB程序；lp=10A(-5);lt=10A(-4);gm=0.01;A=Nlp+lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0;l,v=eig(A);c=inv(1)*1;0;0輸出結果：特征值(0,-0.0102,-0.0001)與之對應的特征向量(001)；(0.7053-0.70890.0035);(-0.7053-0.00350.7089)根據(jù)特征值和特征向量寫出其解析解，利用初始條件田地遍，1，0，0得到特解0.0102t0.0001tP1 (t) = 0.004937

51、e0.99503724eP2(t)-0.0049623e0.0102t00.000t0.0049378eP3(t) =1 0.0000245e00102t -1.00011612e。0001t.數(shù)值解法MATLAB程序如下首先編寫m文件(eqs0. m)function xdot=eqs0(t,p,flag,lp,lt,gm)a=-(lp=lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+ lt,0；p=p(1);p(2);p(3) xdot=a*p;在工作空間執(zhí)行以下程序:ts=01000;p0=1;0;0;lp=10A(-5);lt=10A(-4);gm=0.0

52、1;t,p=ode23(' eqs0' ,ts,p0,lp,lt,gm)plot(t,1-p(:,);xlabel(時間 t(小時);ylabe,靠度 R(t)');400 WO時閶M小mijo .系統(tǒng)正常運行的可靠度曲蛾參數(shù)取偏 gXOOOOLlHOML和尸0title(參數(shù)取值lp=0.00001;lt=0.00001;gm=0.01');grid輸出結果如圖3.10計算表明：當t=1000小時情況下，p，t),p2(t),p3(t)目0.9369,0.0047,0.0584顯示系統(tǒng)工作的概率為94.18%,從圖3.10中還可以看出系統(tǒng)可靠性變化更加細致的情

53、況，何時可靠度達到99%,何時達到98%等等。.敏感性分析H10，S：0.010.1內(nèi)變化時，系統(tǒng)2x2、dxdt2 "oooo-x )&- x(0) =2;x'=0當參數(shù)取值在以下范圍p：1010例14y1 ' = y22、y2' = 1000(1 - y1 )y2 - y1y(0) =2, y2=0的可靠性情況如何呢？留給讀者自己分析解：令yi=x,y2=yi'則微分方程變?yōu)橐浑A微分方程組:1、建立m-文件vdp1000.m如下:functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)A2)*y(2)-y(1);2、取t0=0,tf=3000,輸入命令:T,Y=ode15s('vdp1000',03000,20);plot(T,Y(:,1),'-')3、結果如圖程序見數(shù)學建模高教出版社數(shù)學建模與數(shù)學實驗ToMatlab(ff4)導彈追蹤問題設位于坐標原點的甲艦向位于x軸上點A(1,0)處的乙艦發(fā)射導彈，導彈頭始終對準乙艦.如果乙艦以最大的速度v0(是常數(shù))沿平行于y軸的直線行駛，導彈的速度是5v0,求導彈運行的曲線方程.又乙艦行駛多遠時，導彈將它