微分方程數(shù)值解課程設(shè)計(jì)報(bào)告

1、. -. -可修編- . z.微分方程數(shù)值方法課程設(shè)計(jì)Basic Theory of Ordinary Differential Equations E*periment Report教務(wù)處2014年7月課 程 設(shè) 計(jì) 說(shuō) 明 書題目：常微分方程數(shù)值解法課程設(shè)計(jì)學(xué)院系：理學(xué)院年級(jí)專業(yè)：計(jì)算科學(xué)11-1學(xué)生：指導(dǎo)教師：教師職稱：教授燕山大學(xué)課程設(shè)計(jì)論文任務(wù)書院系：理學(xué)院教學(xué)單位：學(xué) 號(hào)*學(xué)生*專業(yè)班級(jí)11計(jì)算數(shù)學(xué)設(shè)計(jì)題目常微分方程數(shù)值解法課程設(shè)計(jì)設(shè)計(jì)容比擬Adams四階PECE模式和PMECME模式。求解貝塞爾方程并與準(zhǔn)確解比擬。小型火箭初始重量為1400kg，其中包括1080kg 燃料?；鸺Q

2、直向上發(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、要通過(guò)課題背景的介紹，說(shuō)明選題的意義；2、運(yùn)用所學(xué)知識(shí)，建立問(wèn)題的數(shù)學(xué)模型；3、利用Matlab軟件進(jìn)展計(jì)算，并給出計(jì)算機(jī)程序框圖及計(jì)算機(jī)程序；4、運(yùn)用所學(xué)知識(shí)，對(duì)計(jì)算結(jié)果進(jìn)展分析；工作量1、程設(shè)計(jì)要在兩周完成；2、根據(jù)題目找出解題方法，然后利用Matlab進(jìn)展計(jì)算；3、正文在3000字左右。工作計(jì)劃第一天：收集

3、資料，閱讀參考文獻(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é)單位主任簽字說(shuō)明：此表一式四份，學(xué)生、指導(dǎo)教師、基層教學(xué)單位、系部各一份。年 月 日 燕山大學(xué)課程設(shè)計(jì)評(píng)審意見表指導(dǎo)教師評(píng)語(yǔ)：成績(jī)：指導(dǎo)教師：年 月 日辯論小組評(píng)語(yǔ)：成績(jī)：評(píng)閱

4、人：年 月 日課程設(shè)計(jì)總成績(jī)：辯論小組成員簽字：年 月 日-. z.摘 要本文對(duì)常微分方程初值問(wèn)題現(xiàn)有的數(shù)值解法進(jìn)展了綜述研究。主要討論了幾種常用的數(shù)值解法：即歐拉法，改良?xì)W拉法，龍格庫(kù)塔方法，阿達(dá)姆斯PECE,PMECME格式等。文章最后結(jié)合常見數(shù)值解法，對(duì)較為典型的微分方程模型進(jìn)展數(shù)值求解，探討了上述數(shù)值算法在實(shí)際建模問(wèn)題中的應(yīng)用。論文闡述的是常微分方程數(shù)值解法的幾個(gè)問(wèn)題，通過(guò)對(duì)以下問(wèn)題的求解一.比擬Adams四階PECE模式和PMECME模式。二.求解貝塞爾方程并與準(zhǔn)確解比擬。三.小型火箭初始重量為1400kg，其中包括1080kg 燃料?；鸺Q直向上發(fā)射時(shí)燃料燃燒率為18kg/s，由此

5、產(chǎn)生32000N 的推力，火箭引擎在燃料用盡時(shí)關(guān)閉。設(shè)火箭上升時(shí)空氣阻力正比于速度的平方，比例系數(shù)為0.4kg/m，求引擎關(guān)閉瞬間火箭的高度、速度、加速度，及火箭到達(dá)最高點(diǎn)時(shí)的高度和加速度，并畫出高度、速度、加速度隨時(shí)間變化的圖形。來(lái)加強(qiáng)對(duì)用數(shù)值解法解常微分方程實(shí)際問(wèn)題的能力。關(guān)鍵詞：常微分方程數(shù)值解 MATLAB-. z.目 錄 TOC o 1-3 h z HYPERLINK l _Toc293560829摘要 PAGEREF _Toc293560829 h iHYPERLINK l _Toc293560830緒論 PAGEREF _Toc293560830 h 1問(wèn)題解答.2HYPERLI

6、NK l _Toc293560842總結(jié).17HYPERLINK l _Toc293560843參考文獻(xiàn)資料.17. -. z.緒 論很多科學(xué)技術(shù)和工程問(wèn)題常用微分方程的形式建立數(shù)學(xué)模型，因此微分方程的求解是很有意義的。建立微分方程只是解決問(wèn)題的第一步，通常需要求出方程的解來(lái)說(shuō)明實(shí)際現(xiàn)象，并加以檢驗(yàn)。如果能得到解析形式的解固然是便于分析和應(yīng)用的，但是我們知道，雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來(lái)求解一些典型的方程，而對(duì)于絕大多數(shù)的微分方程問(wèn)題，很難或者根本不可能得到它的解析解，實(shí)際問(wèn)題終歸結(jié)出來(lái)的微分方程主要靠數(shù)值解法。因此，研究微分方程求解的數(shù)值方法是非常有意義的。-.

7、 z.問(wèn)題解答問(wèn)題一：比擬Adams四階PECE模式和PMECME模式。問(wèn)題引出：將阿達(dá)姆斯方法顯式與隱式方法作一比照，以說(shuō)明預(yù)測(cè)校正格式的必要性。這些方法的階及誤差常數(shù)列表如下：步數(shù)階誤差常數(shù)顯式隱式顯式隱式112223334445由于阿達(dá)姆斯插公式是隱式公式，故用它計(jì)算時(shí)也需用迭代法。通常把阿達(dá)姆斯外插公式與插公式結(jié)合起來(lái)使用，先由前者提供初值，再由后者進(jìn)展修正，所以Adams預(yù)測(cè)-校正格式既利用了隱式方法較好的穩(wěn)定性及準(zhǔn)確性，有利用了顯式方法的簡(jiǎn)易性，正是把兩者結(jié)合起來(lái)，做到取長(zhǎng)補(bǔ)短。Adams四階預(yù)估-校正(PECE)公式： 而PMECME模式公式為對(duì)于初值問(wèn)題u=u-2t/u,u(0

8、)=1,在單區(qū)間【0,1】上，用Adams四階預(yù)估-校正算法的PECE模式及PMECME模式求其數(shù)值解，取步長(zhǎng)h=0.1利用計(jì)算結(jié)果估計(jì)數(shù)值解的局部誤差主項(xiàng)。真解ut=sqrt(1+2t)編寫matlab程序：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=

9、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(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

10、*T(k)/U(k);endUn=U;for k=1:11 zz(k)=sqrt(1+2*T(k);ende=U-zz;endfunction Un,e = pmecme()% Un,e = pmecme()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

11、/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(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)

12、-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é)果及圖形顯示：運(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

13、 0.0571 0.0271 0.0127 0.0042 Columns 9 through 11 -0.0013 -0.0054 -0.0088pmecmeans = Columns 1 through 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)

14、plot(t1,U1,-r,Marker,*)plot(t1,U2,-g,Marker,*)legend(U=sqrt(1+2*t),U1,U2Location,NorthWest)grid on*lable(t)ylable(Un)hold off運(yùn)行結(jié)果如下：放大如下列圖可知PMECME模式比PECE模式更為準(zhǔn)確問(wèn)題二：求解貝塞爾方程并與準(zhǔn)確解比擬。求解 *2y+*y+(*2-n2)y=0 y=2,y=- (貝賽爾方程，令n=0.5)，準(zhǔn)確解y=sin *解：首先將高階方程裝降階化簡(jiǎn)為一階常微分方程組令y1=y，令y2=y1=y將n=0.5代入，則原方程轉(zhuǎn)化為：y1=y2y1=2, y2=

15、-1用向前歐拉公式：y1(n+1)=y1(n)+h*y2(n)y2(n+1)=y2(n)+h*y1(0) =2，y2(0)= -，*(n+1)=+n*h，h為步長(zhǎng)編程如下：clear all*=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)/*(i)+(1-0.25/*(i)2)*y1(i);endplot(*,y1),grid但如果步長(zhǎng)選得不一樣(橫坐標(biāo)都從開場(chǎng)，到100左右完畢，但中間選的點(diǎn)也不太一樣)，即使采用同樣的迭代公式，繪出的曲線卻有

16、很大不同：程序：clear all*=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(i+1)=y2(i)-0.01*pi*(y2(i)/*(i)+(1-0.25/*(i)2)*y1(i);endplot(*,y1),grid我們認(rèn)為，可能是因?yàn)闅W拉公式每算一次都會(huì)產(chǎn)生誤差，如果取的點(diǎn)正好位置不太適宜，會(huì)導(dǎo)致誤差一步步增大，累加起來(lái)，就像蝴蝶效應(yīng)一樣，會(huì)產(chǎn)生和真實(shí)狀況截然不同的結(jié)果。這也是用數(shù)值方法求解方程最怕出現(xiàn)的問(wèn)題，也是應(yīng)該解決的問(wèn)題。這說(shuō)明向前歐拉公式并不是一種很好的計(jì)算方法，誤

17、差較大只有1階精度，而且容易失真。這一點(diǎn)在下面還要說(shuō)明。2龍格庫(kù)塔方法y1=y2y1=2, y2=-編寫如下M文件：function d*=fangchengzu(t,*) %建立名為fangchengzu的函數(shù)M文件d*=*(2);-*(2)/t-(1-0.25/t2)*(1); %表示方程組用4級(jí)5階龍格庫(kù)塔公式進(jìn)展計(jì)算。編程如下：clear allts=pi/2:0.1:100;%使用龍格庫(kù)塔方法*0=2,-2/pi;t,*=ode45(fangchengzu,ts,*0);plot(t,*(:,1),grid,title(龍格-庫(kù)塔方法),pause%繪出準(zhǔn)確解y=sin *圖像t=p

18、i/2:0.1:100;y=sin(t).*sqrt(2*pi./t);plot(t,y),grid,title(準(zhǔn)確解)繪出曲線如下：比擬分析：從曲線上可以看出，在數(shù)值求解貝賽爾方程時(shí)，龍格庫(kù)塔方法4級(jí)5階龍格庫(kù)塔公式的結(jié)果與準(zhǔn)確解y=sin *非常接近，產(chǎn)生震蕩，并且隨*的增大振幅逐漸減?。欢脷W拉方法向前歐拉公式的計(jì)算結(jié)果卻與準(zhǔn)確解有較大差異，雖然也產(chǎn)生了震蕩，但在*稍大的情況下，振幅隨*的增大而增大。這是因?yàn)橄蚯皻W拉公式只有1階的精度，而4級(jí)5階龍格庫(kù)塔公式有5階的精度，一般情況下，肯定用后者計(jì)算更接近準(zhǔn)確解。而且從本例可以看出，向前歐拉公式并不是一種很好的計(jì)算方法，誤差較大，甚至導(dǎo)致

19、結(jié)果失真，所以還是應(yīng)該用龍格庫(kù)塔方法比擬好。問(wèn)題三：小型火箭初始質(zhì)量為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í)間變化的圖形。1、簡(jiǎn)要分析此題的求解需要用到常微分方程，而整個(gè)過(guò)程又被分為兩個(gè)階段：火箭加速上升階段和燃料燃盡后減速的階段。由題目易知第一個(gè)階段持續(xù)時(shí)間T1=60s列出第一階段的方程組：設(shè)M0為火箭本身質(zhì)量，m為燃料質(zhì)量，g為重

20、力加速度 = 9.8m/,燃料燃燒率為a，空氣阻力的比例系數(shù)為k，F(xiàn)為推進(jìn)力。M0 = 1400-1080 = 320kg; V=(F-kv2-M0+mg)/(M0+m) QUOTE m=-a QUOTE 初值 QUOTE v=0,m=1080。由以上各式可以求出t= QUOTE T1時(shí)火箭的速度。再求解第二階段： V=(-kv2-M0g)/M0 QUOTE m=0 QUOTE 可以求出火箭速度降為0的時(shí)刻。將整個(gè)過(guò)程中的時(shí)間向量以及速度向量聯(lián)合起來(lái)，利用第三章所學(xué)插值與數(shù)值積分的方法可以求得任意時(shí)刻火箭的近似高度。2、方法求解常微分方程時(shí)，我分別采用了自己編寫的歐拉公式、改良?xì)W拉公式、4級(jí)4

21、階龍格-庫(kù)塔公式，以及MATLAB自帶的龍格-庫(kù)塔方法，結(jié)果與分析1、各種公式的比照首先，我作出了各種不同公式計(jì)算得到的火箭速度隨時(shí)間變化的圖像，圖如下：從圖中可以看出，各種公式計(jì)算得到的結(jié)果根本一致，為確定其區(qū)別，將圖像放大，放大約2000 倍后，得到下列圖：分析：從圖中可以看出，自編歐拉公式距離MATLAB自帶龍格-庫(kù)塔公式最遠(yuǎn)，精度最差；自編的改良?xì)W拉公式和自編的龍格-庫(kù)塔公式結(jié)果根本一致，兩者中自編龍格-庫(kù)塔公式距MATLAB自帶龍格-庫(kù)塔公式的結(jié)果稍近。與之前的分析根本一致。然而產(chǎn)生自編龍格-庫(kù)塔公式與MATLAB自帶龍格-庫(kù)塔公式之間的差距的原因還未知。由于MATLAB自帶龍格-庫(kù)

22、塔公式精度較高，因此以下不在計(jì)算其它幾項(xiàng)公式的結(jié)果。2、第一階段火箭關(guān)閉瞬間的速度 QUOTE v=267.1m/s關(guān)閉前瞬間的加速度a=(F -kv2-M0g)/M0=0.1320m/s2 QUOTE 此時(shí)火箭的高度h=12189.7791507305m3、第二階段由初始速度以及常微分方程可以求得火箭到達(dá)最高點(diǎn)的時(shí)間約為71.3s；此時(shí)火箭的高度 QUOTE h=13115.748739887m加速度a=-9.111m/s24、整體過(guò)程下列圖為火箭加速度與時(shí)間的關(guān)系。由圖可以看出，火箭一開場(chǎng)的加速度很大，隨著時(shí)間的推移，火箭的燃料有所減少，與此同時(shí)速度有所上升。兩者中前者使火箭的加速度增大，

23、后者使其減小，綜合作用，最終表達(dá)為加速度先略有上升，然后慢慢減小。當(dāng)火箭中燃料燃盡時(shí)，火箭喪失推動(dòng)力，因而加速度急劇減小為負(fù)值。此后火箭速度不斷減小，導(dǎo)致火箭所受阻力逐漸減小，因而加速度的絕對(duì)值有所減小，直到最終火箭速度降為零，火箭不受阻力，僅受重力，加速度為重力加速度。下列圖為火箭速度與時(shí)間的關(guān)系：此圖可以看作是由第一幅圖對(duì)時(shí)間積分所得結(jié)果，本圖的斜率對(duì)應(yīng)第一幅圖的值。第一階段火箭加速度為正，因此速度不斷增加，只是增加的速度不斷減慢。燃料燃盡后，加速度變?yōu)樨?fù)值，因此速度開場(chǎng)急劇下降，與此同時(shí)下降的速率不斷減小。最終火箭速度降為0。下列圖為火箭高度與時(shí)間的關(guān)系：此圖可以看作是由第二幅圖對(duì)時(shí)間積

24、分所得結(jié)果，本圖的斜率對(duì)應(yīng)第一幅圖的值。一開場(chǎng)或減速度較小，因此高度緩慢增加，之后增加的速度不斷提升，直到火箭的燃料耗盡，此后速度不斷減小，但仍為正值，因此火箭繼續(xù)向上飛行，只是高度提升的速度逐漸變慢。知道最終火箭速度降為零。程序清單1、第一階段的微分方程組function d* = fly(t,y)M0 = 320;a = 18; k = 0.4;g = 9.8;F = 32000;d* = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; - a*sign(y(2);2、第二階段微分方程組function d* = fly1(t,y)M0 = 320; k

25、= 0.4;g = 9.8;F = 0;a = 0;d* = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; -a;3、自編歐拉公式function y = Euler(fun,ts,*0)NumOfLine = numel(ts);NumOfRow = numel(*0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = *0;for i = 1:NumOfLine-1 dy = feval(fun,ts(i),y(i,:); y(i+1,:)=y(i,:)+rot90(dy)*h;end4、自編改

26、良?xì)W拉公式function y = ImprovedEuler(fun,ts,*0)NumOfLine = numel(ts);NumOfRow = numel(*0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = *0;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、自編龍格-庫(kù)塔公式fun

27、ction y = myRungeKutta(fun,ts,*0)NumOfLine = numel(ts);NumOfRow = numel(*0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = *0;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

28、); 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 =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),

29、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