微分方程數(shù)值解課程設(shè)計(jì)

1、 微分方程數(shù)值方法課程設(shè)計(jì)Basic Theory of Ordinary Differential Equations Experiment Report教 務(wù) 處2014年 7 月課 程 設(shè) 計(jì) 說 明 書題目： 常微分方程數(shù)值解法課程設(shè)計(jì) 學(xué)院（系）： 理學(xué)院 年級專業(yè)： 計(jì)算科學(xué)11-1 學(xué)生姓名： 指導(dǎo)教師： 教師職稱： 教授 燕山大學(xué)課程設(shè)計(jì)（論文）任務(wù)書院（系）：理學(xué)院 教學(xué)單位： 學(xué) 號XXX學(xué)生姓名XXX專業(yè)（班級）11計(jì)算數(shù)學(xué)設(shè)計(jì)題目常微分方程數(shù)值解法課程設(shè)計(jì)設(shè)計(jì)內(nèi)容一. 比較Adams四階PECE模式和PMECME模式。二. 求解貝塞爾方程并與精確解比較。三. 小型火箭初

2、始重量為1400kg，其中包括1080kg 燃料?；鸺Q直向上發(fā)射時(shí)燃料燃燒率為18kg/s，由此產(chǎn)生32000N 的推力，火箭引擎在燃料用盡時(shí)關(guān)閉。設(shè)火箭上升時(shí)空氣阻力正比于速度的平方，比例系數(shù)為0.4kg/m，求引擎關(guān)閉瞬間火箭的高度、速度、加速度，及火箭到達(dá)最高點(diǎn)時(shí)的高度和加速度，并畫出高度、速度、加速度隨時(shí)間變化的圖形。設(shè)計(jì)要求1、要通過課題背景的介紹，闡明選題的意義；2、運(yùn)用所學(xué)知識，建立問題的數(shù)學(xué)模型；3、利用Matlab軟件進(jìn)行計(jì)算，并給出計(jì)算機(jī)程序框圖及計(jì)算機(jī)程序；4、運(yùn)用所學(xué)知識，對計(jì)算結(jié)果進(jìn)行分析；工作量1、程設(shè)計(jì)要在兩周內(nèi)完成；2、根據(jù)已知題目找出解題方法，然后利用Mat

3、lab進(jìn)行計(jì)算；3、正文在3000字左右。工作計(jì)劃第一天：收集資料，閱讀參考文獻(xiàn)；第二天：構(gòu)建數(shù)學(xué)模型；第三天：設(shè)計(jì)計(jì)算機(jī)流程圖，編制計(jì)算機(jī)程序，上機(jī)運(yùn)算，調(diào)試程序；第四天：計(jì)算結(jié)果分析第五天：撰寫報(bào)告等參考資料1胡建偉湯懷民。微分方程數(shù)值方法?？茖W(xué)出版社。1999年01月2 陳渝周璐錢方等。數(shù)值方法（MATLAB版）。電子工業(yè)出版社。2002.63尹澤明，丁春麗等。精通MATLAB 6。清華大學(xué)出版社。2006.64Cleve B.Moler。MATLAB數(shù)值計(jì)算。機(jī)械工業(yè)出版社2006.6指導(dǎo)教師簽字基層教學(xué)單位主任簽字說明：此表一式四份，學(xué)生、指導(dǎo)教師、基層教學(xué)單位、系部各一份。年 月

4、日 燕山大學(xué)課程設(shè)計(jì)評審意見表指導(dǎo)教師評語：成績： 指導(dǎo)教師： 年 月 日答辯小組評語：成績： 評閱人： 年 月 日課程設(shè)計(jì)總成績：答辯小組成員簽字：年 月 日摘 要本文對常微分方程初值問題現(xiàn)有的數(shù)值解法進(jìn)行了綜述研究。主要討論了幾種常用的數(shù)值解法：即歐拉法，改進(jìn)歐拉法，龍格庫塔方法，阿達(dá)姆斯PECE,PMECME格式等。文章最后結(jié)合常見數(shù)值解法，對較為典型的微分方程模型進(jìn)行數(shù)值求解，探討了上述數(shù)值算法在實(shí)際建模問題中的應(yīng)用。論文闡述的是常微分方程數(shù)值解法的幾個(gè)問題，通過對以下問題的求解一.比較Adams四階PECE模式和PMECME模式。二.求解貝塞爾方程并與精確解比較。三.小型火箭初始重量

5、為1400kg，其中包括1080kg 燃料。火箭豎直向上發(fā)射時(shí)燃料燃燒率為18kg/s，由此產(chǎn)生32000N 的推力，火箭引擎在燃料用盡時(shí)關(guān)閉。設(shè)火箭上升時(shí)空氣阻力正比于速度的平方，比例系數(shù)為0.4kg/m，求引擎關(guān)閉瞬間火箭的高度、速度、加速度，及火箭到達(dá)最高點(diǎn)時(shí)的高度和加速度，并畫出高度、速度、加速度隨時(shí)間變化的圖形。來加強(qiáng)對用數(shù)值解法解常微分方程實(shí)際問題的能力。關(guān)鍵詞：常微分方程數(shù)值解 MATLAB17目 錄摘 要i緒 論1問題解答.2總結(jié).17參考文獻(xiàn)資料.171 緒論緒 論很多科學(xué)技術(shù)和工程問題常用微分方程的形式建立數(shù)學(xué)模型，因此微分方程的求解是很有意義的。建立微分方程只是解決問題的

6、第一步，通常需要求出方程的解來說明實(shí)際現(xiàn)象，并加以檢驗(yàn)。如果能得到解析形式的解固然是便于分析和應(yīng)用的，但是我們知道，雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些典型的方程，而對于絕大多數(shù)的微分方程問題，很難或者根本不可能得到它的解析解，實(shí)際問題終歸結(jié)出來的微分方程主要靠數(shù)值解法。因此，研究微分方程求解的數(shù)值方法是非常有意義的。燕 山 大 學(xué) 課 程 設(shè) 計(jì) 說 明 書問題解答問題一：比較Adams四階PECE模式和PMECME模式。問題引出：將阿達(dá)姆斯方法顯式與隱式方法作一對比，以說明預(yù)測校正格式的必要性。這些方法的階及誤差常數(shù)列表如下：步數(shù)階誤差常數(shù)顯式隱式顯式隱式11

7、2223334445由于阿達(dá)姆斯內(nèi)插公式是隱式公式，故用它計(jì)算時(shí)也需用迭代法。通常把阿達(dá)姆斯外插公式與內(nèi)插公式結(jié)合起來使用，先由前者提供初值，再由后者進(jìn)行修正，所以Adams預(yù)測-校正格式既利用了隱式方法較好的穩(wěn)定性及精確性，有利用了顯式方法的簡易性，正是把兩者結(jié)合起來，做到取長補(bǔ)短。Adams四階預(yù)估-校正(PECE)公式： 而PMECME模式公式為對于初值問題u=u-2t/u,u(0)=1,在單區(qū)間【0,1】上，用Adams四階預(yù)估-校正算法的PECE模式及PMECME模式求其數(shù)值解，取步長h=0.1利用計(jì)算結(jié)果估計(jì)數(shù)值解的局部誤差主項(xiàng)。（真解u（t）=sqrt(1+2t)）編寫matla

8、b程序：function Un,e = pece()% Un,e = pece()f0=u-2*t/u;v=t,u;U=zeros(1,11);T=0:0.1:1;f=zeros(1,11);h=0.1;U(1)=1;f(1)=1;for k=1:3 K1=subs(f0,v,(k-1)*h,U(k); K2=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K1); K3=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K2); K4=subs(f0,v,(k-1)*h+h,U(k)+h*K3); U(k+1)=U(k)+1/6*h*(K1+2*K2

9、+2*K3+K4); f(k+1)=U(k+1)-2*T(k+1)/U(k+1);endfor k=5:11 U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4); f(k)=U(k)-2*T(k)/U(k); U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3); f(k)=U(k)-2*T(k)/U(k);endUn=U;for k=1:11 zz(k)=sqrt(1+2*T(k);ende=U-zz;end function Un,e = pmecme()% Un,e = pmec

10、me()syms t uf0=u-2*t/u;v=t,u;U=zeros(1,11);T=0:0.1:1;f=zeros(1,11);h=0.1;U(1)=1;f(1)=1;for k=1:3 K1=subs(f0,v,(k-1)*h,U(k); K2=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K1); K3=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K2); K4=subs(f0,v,(k-1)*h+h,U(k)+h*K3); U(k+1)=U(k)+1/6*h*(K1+2*K2+2*K3+K4); f(k+1)=U(k+1)-2*T(

11、k+1)/U(k+1);endfor k=5:11 U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4); U(k)=U(k)+251/270*(U(k)-U(k-1); f(k)=U(k)-2*T(k)/U(k); U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3); U(k)=U(k)-19/270*(U(k)-U(k-1); f(k)=U(k)-2*T(k)/U(k);endUn=U;for k=1:11 zz(k)=sqrt(1+2*T(k);ende=U-zz;end 運(yùn)行結(jié)

12、果及圖形顯示：運(yùn)行得： Un,e=pece()Un = Columns 1 through 8 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 Columns 9 through 11 1.6125 1.6733 1.7321e = 1.0e-05 * Columns 1 through 8 0 0.0417 0.0789 0.1164 0.0571 0.0271 0.0127 0.0042 Columns 9 through 11 -0.0013 -0.0054 -0.0088pmecmeans = Columns 1 thro

13、ugh 8 1.0000 1.0954 1.1832 1.2649 1.3398 1.4104 1.4773 1.5409 Columns 9 through 11 1.6016 1.6595 1.7147編寫a.m用于畫圖顯示計(jì)算結(jié)果，t0=0:0.001:1;t1=0:0.1:1;U0=sqrt(1+2*t0);U2=pece();U1=pmecme();hold onplot(t0,U0,k)plot(t1,U1,-r,Marker,x)plot(t1,U2,-g,Marker,x)legend(U=sqrt(1+2*t),U1,U2Location,NorthWest)grid onx

14、lable(t)ylable(Un)hold off運(yùn)行結(jié)果如下：放大如下圖可知PMECME模式比PECE模式更為精確問題二：求解貝塞爾方程并與精確解比較。求解 x2y+xy+(x2-n2)y=0 y=2,y=- (貝賽爾方程，令n=0.5)，精確解y=sin x解：首先將高階方程裝降階化簡為一階常微分方程組令y1=y，令y2=y1=y將n=0.5代入，則原方程轉(zhuǎn)化為：y1=y2y1=2, y2=-（1）用向前歐拉公式：y1(n+1)=y1(n)+h*y2(n)y2(n+1)=y2(n)+h*y1(0) =2，y2(0)= -，x(n+1)= +n*h，h為步長編程如下：clear allx=

15、pi/2:0.1:pi/2+100-0.05;y1=2;y2=-2/pi;for i=1:999 y1(i+1)=y1(i)+0.1*y2(i); y2(i+1)=y2(i)-0.1*(y2(i)/x(i)+(1-0.25/x(i)2)*y1(i);endplot(x,y1),grid但如果步長選得不一樣(橫坐標(biāo)都從開始，到100左右結(jié)束，但中間選的點(diǎn)也不太一樣)，即使采用同樣的迭代公式，繪出的曲線卻有很大不同：程序：clear allx=pi/2:0.01*pi:30*pi;y1=2;y2=-2/pi;for i=1:2950 y1(i+1)=y1(i)+0.01*pi*y2(i); y2(

16、i+1)=y2(i)-0.01*pi*(y2(i)/x(i)+(1-0.25/x(i)2)*y1(i);endplot(x,y1),grid我們認(rèn)為，可能是因?yàn)闅W拉公式每算一次都會(huì)產(chǎn)生誤差，如果取的點(diǎn)正好位置不太合適，會(huì)導(dǎo)致誤差一步步增大，累加起來，就像蝴蝶效應(yīng)一樣，會(huì)產(chǎn)生和真實(shí)狀況截然不同的結(jié)果。這也是用數(shù)值方法求解方程最怕出現(xiàn)的問題，也是應(yīng)該解決的問題。這說明向前歐拉公式并不是一種很好的計(jì)算方法，誤差較大（只有1階精度），而且容易失真。這一點(diǎn)在下面還要說明。（2）龍格庫塔方法y1=y2y1=2, y2=-編寫如下M文件：function dx=fangchengzu(t,x) %建立名為f

17、angchengzu的函數(shù)M文件dx=x(2);-x(2)/t-(1-0.25/t2)*x(1); %表示方程組用4級5階龍格庫塔公式進(jìn)行計(jì)算。編程如下：clear allts=pi/2:0.1:100;%使用龍格庫塔方法x0=2,-2/pi;t,x=ode45(fangchengzu,ts,x0);plot(t,x(:,1),grid,title(龍格-庫塔方法),pause%繪出精確解y=sin x圖像t=pi/2:0.1:100;y=sin(t).*sqrt(2*pi./t);plot(t,y),grid,title(精確解)繪出曲線如下：比較分析：從曲線上可以看出，在數(shù)值求解貝賽爾方程

18、時(shí)，龍格庫塔方法（4級5階龍格庫塔公式）的結(jié)果與精確解y=sin x非常接近，產(chǎn)生震蕩，并且隨x的增大振幅逐漸減小；而用歐拉方法（向前歐拉公式）的計(jì)算結(jié)果卻與精確解有較大差異，雖然也產(chǎn)生了震蕩，但在x稍大的情況下，振幅隨x的增大而增大。這是因?yàn)橄蚯皻W拉公式只有1階的精度，而4級5階龍格庫塔公式有5階的精度，一般情況下，肯定用后者計(jì)算更接近精確解。而且從本例可以看出，向前歐拉公式并不是一種很好的計(jì)算方法，誤差較大，甚至導(dǎo)致結(jié)果失真，所以還是應(yīng)該用龍格庫塔方法比較好。問題三：小型火箭初始質(zhì)量為1400kg，其中包括1080kg燃料，火箭豎直向上發(fā)射時(shí)燃料燃燒率為18kg/s，由此產(chǎn)生32000N的

19、推力，火箭引擎在燃料用盡時(shí)關(guān)閉，設(shè)火箭上升時(shí)空氣阻力正比于速度的平方，比例系數(shù)為0.4kg/m，求引擎關(guān)閉瞬間火箭的高度，速度，加速度及火箭達(dá)到最高點(diǎn)時(shí)的高度和加速度，并畫出高度、速度、加速度隨時(shí)間變化的圖形。1、簡要分析本題的求解需要用到常微分方程，而整個(gè)過程又被分為兩個(gè)階段：火箭加速上升階段和燃料燃盡后減速的階段。由題目易知第一個(gè)階段持續(xù)時(shí)間T1=60s列出第一階段的方程組：設(shè)M0為火箭本身質(zhì)量，m為燃料質(zhì)量，g為重力加速度 = 9.8m/s2,燃料燃燒率為a，空氣阻力的比例系數(shù)為k，F(xiàn)為推進(jìn)力。M0 = 1400-1080 = 320kg; V=(F-kv2-M0+mg)/(M0+m)

20、m=-a初值v=0,m=1080。由以上各式可以求出t=T1時(shí)火箭的速度。再求解第二階段： V=(-kv2-M0g)/M0 m=0可以求出火箭速度降為0的時(shí)刻。將整個(gè)過程中的時(shí)間向量以及速度向量聯(lián)合起來，利用第三章所學(xué)插值與數(shù)值積分的方法可以求得任意時(shí)刻火箭的近似高度。2、方法求解常微分方程時(shí)，我分別采用了自己編寫的歐拉公式、改進(jìn)歐拉公式、4級4階龍格-庫塔公式，以及MATLAB自帶的龍格-庫塔方法，結(jié)果與分析1、各種公式的對比首先，我作出了各種不同公式計(jì)算得到的火箭速度隨時(shí)間變化的圖像，圖如下：從圖中可以看出，各種公式計(jì)算得到的結(jié)果基本一致，為確定其區(qū)別，將圖像放大，放大約2000 倍后，得

21、到下圖：分析：從圖中可以看出，自編歐拉公式距離MATLAB自帶龍格-庫塔公式最遠(yuǎn)，精度最差；自編的改進(jìn)歐拉公式和自編的龍格-庫塔公式結(jié)果基本一致，兩者中自編龍格-庫塔公式距MATLAB自帶龍格-庫塔公式的結(jié)果稍近。與之前的分析基本一致。然而產(chǎn)生自編龍格-庫塔公式與MATLAB自帶龍格-庫塔公式之間的差距的原因還未知。由于MATLAB自帶龍格-庫塔公式精度較高，因此以下不在計(jì)算其它幾項(xiàng)公式的結(jié)果。2、第一階段火箭關(guān)閉瞬間的速度v=267.261240773261m/s關(guān)閉前瞬間的加速度a=(F -kv2-M0g)/M0=0.914286475421320m/s2此時(shí)火箭的高度h=12189.77

22、91507305m3、第二階段由初始速度以及常微分方程可以求得火箭達(dá)到最高點(diǎn)的時(shí)間約為71.3s；此時(shí)火箭的高度h=13115.748739887m加速度a=-9.80000406052111m/s24、整體過程下圖為火箭加速度與時(shí)間的關(guān)系。由圖可以看出，火箭一開始的加速度很大，隨著時(shí)間的推移，火箭的燃料有所減少，與此同時(shí)速度有所上升。兩者中前者使火箭的加速度增大，后者使其減小，綜合作用，最終體現(xiàn)為加速度先略有上升，然后慢慢減小。當(dāng)火箭中燃料燃盡時(shí)，火箭喪失推動(dòng)力，因而加速度急劇減小為負(fù)值。此后火箭速度不斷減小，導(dǎo)致火箭所受阻力逐漸減小，因而加速度的絕對值有所減小，直到最終火箭速度降為零，火箭

23、不受阻力，僅受重力，加速度為重力加速度。 下圖為火箭速度與時(shí)間的關(guān)系： 此圖可以看作是由第一幅圖對時(shí)間積分所得結(jié)果，本圖的斜率對應(yīng)第一幅圖的值。第一階段火箭加速度為正，因此速度不斷增加，只是增加的速度不斷減慢。燃料燃盡后，加速度變?yōu)樨?fù)值，因此速度開始急劇下降，與此同時(shí)下降的速率不斷減小。最終火箭速度降為0。 下圖為火箭高度與時(shí)間的關(guān)系：此圖可以看作是由第二幅圖對時(shí)間積分所得結(jié)果，本圖的斜率對應(yīng)第一幅圖的值。一開始或減速度較小，因此高度緩慢增加，之后增加的速度不斷提升，直到火箭的燃料耗盡，此后速度不斷減小，但仍為正值，因此火箭繼續(xù)向上飛行，只是高度提升的速度逐漸變慢。知道最終火箭速度降為零。程序

24、清單1、第一階段的微分方程組function dx = fly(t,y)M0 = 320;a = 18; k = 0.4;g = 9.8;F = 32000;dx = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; - a*sign(y(2);2、第二階段微分方程組function dx = fly1(t,y)M0 = 320; k = 0.4;g = 9.8;F = 0;a = 0;dx = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; -a;3、自編歐拉公式function y = Euler(fun,ts,x0)Num

25、OfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = x0;for i = 1:NumOfLine-1 dy = feval(fun,ts(i),y(i,:); y(i+1,:)=y(i,:)+rot90(dy)*h;end4、自編改進(jìn)歐拉公式function y = ImprovedEuler(fun,ts,x0)NumOfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumO

26、fLine,NumOfRow);y(1,:) = x0;for i = 1:NumOfLine-1 dy = feval(fun,ts(i),y(i,:); y1 = y(i,:)+rot90(dy)*h; dy1 = feval(fun,ts(i)+h,y1); y(i+1,:)=y(i,:)+h/2*(rot90(dy+dy1);end5、自編龍格-庫塔公式function y = myRungeKutta(fun,ts,x0)NumOfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumOfLine,Num

27、OfRow);y(1,:) = x0;for i = 1:NumOfLine-1 k1 = feval(fun,ts(i),y(i,:); dy1 = y(i,:)+h/2*rot90(k1); k2 = feval(fun,ts(i)+h/2,dy1); dy2 = y(i,:)+h/2*rot90(k2); k3 = feval(fun,ts(i)+h/2,dy2); dy3 = y(i,:)+h*rot90(k3); k4 = feval(fun,ts(i)+h,dy3); y(i+1,:)=y(i,:)+h/6*rot90(k1+2*k2+2*k3+k4);end6、“自啟動(dòng)”腳本1e

28、 =0.1;ts1 = 0 : e : 60;y00 = 0,1080;y01 = Euler(fly,ts1,y00);y11 = ImprovedEuler(fly,ts1,y00);y21 = myRungeKutta(fly,ts1,y00);t1,y1 = ode45(fly,ts1,y00);ts2 = 60 : e :80;y00 = y1(60/e+1),1),0;y000 = y01(60/e+1),1),0;y001 = y11(60/e+1),1),0;y002 = y21(60/e+1),1),0;y02 = Euler(fly1,ts2,y000);y12 = ImprovedEuler(fly1,ts1,y001);y22 = myRungeKutta(fly1,ts1,y002);t2,y2 = ode45(fly1,ts2,y00);ts = ts1,ts2; t = t1;t2;i=1;while(y2(i,1)0) i = i+1;endy = y1;y2(1:i-1,:);y0 = y01;y02(1:i-1,:);y1 = y11;y12(1:i-1,:);y2 = y21;y22(1:i-1,:);j = numel(y)/numel(y00);plot(t(1:j),y(:,1),k,ts(1:j),y0(:,1),b,ts(