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

1、<p><b>　　數(shù)值分析課程設(shè)計(jì)</b></p><p>　　題目：常微分方程的數(shù)值解法</p><p><b>　　一、摘要3</b></p><p><b>　　二、設(shè)計(jì)目的3</b></p><p><b>　　三、理論基礎(chǔ)3</b&g

2、t;</p><p>　　1.歐拉公式【3】3</p><p>　　2.改進(jìn)Euler方法【3】3</p><p>　　3.三階龍格-庫(kù)塔方法【3】3</p><p>　　4.四階龍格-庫(kù)塔方法【3】4</p><p>　　四、程序代碼及運(yùn)算結(jié)果4</p><p>　　1、用歐拉法求解

3、5</p><p>　　2、用改進(jìn)歐拉法求解：6</p><p>　　3、用3階龍格—庫(kù)塔求解7</p><p>　　4、 用4階龍格—庫(kù)塔求解9</p><p>　　6、歐拉方法與改進(jìn)歐拉方法、3階龍格－庫(kù)塔法以及4階龍格－庫(kù)塔法得出的解得比較。10</p><p>　　五、數(shù)值分析設(shè)計(jì)的gui界面13&

4、lt;/p><p><b>　　六、結(jié)果分析14</b></p><p><b>　　七、設(shè)計(jì)心得14</b></p><p><b>　　八、參考文獻(xiàn)14</b></p><p><b>　　一、摘要</b></p><p>　

5、　在matlab環(huán)境下熟悉的運(yùn)用計(jì)算機(jī)編程語(yǔ)言并結(jié)合龍格－庫(kù)塔法的理論基礎(chǔ)對(duì)常微分方程初值問(wèn)題進(jìn)行求解，在運(yùn)行完程序后以及對(duì)運(yùn)行結(jié)果做出各方面的分析和比較。</p><p><b>　　二、設(shè)計(jì)目的</b></p><p>　　用熟悉的計(jì)算機(jī)語(yǔ)言編程上機(jī)完成用歐拉方法、改進(jìn)歐拉方法、3階龍格－庫(kù)塔法以及4階龍格－庫(kù)塔法求解常微分方程初值問(wèn)題。</p>&

6、lt;p><b>　　三、理論基礎(chǔ)</b></p><p><b>　　1.歐拉公式【3】</b></p><p>　　在點(diǎn) 將作Taylor展開(kāi)，得</p><p>　　，那么當(dāng)h充分小時(shí)，略去誤差項(xiàng)，用近似替代、近似替代，并注意到，便得</p><p>　　上述方法稱(chēng)為Euler方法。&l

7、t;/p><p>　　2.改進(jìn)Euler方法【3】</p><p>　　在應(yīng)用梯形方法的迭代公式進(jìn)行運(yùn)算時(shí)，每迭代一次都要重新計(jì)算函數(shù)的值，且還要判斷何時(shí)可以終止或轉(zhuǎn)下一步計(jì)算。為了控制計(jì)算量和簡(jiǎn)化算法，通常只迭代一次就轉(zhuǎn)入下一步計(jì)算。具體說(shuō)，我們先用Euler公式求得一個(gè)初步的近似值，稱(chēng)之為預(yù)測(cè)值，然后用梯形方法的迭代公式作一次迭代得，即將校正一次，這樣建立的預(yù)測(cè)-校正方法稱(chēng)之為改進(jìn)的Eul

8、er方法：</p><p><b>　　預(yù)測(cè)：，</b></p><p><b>　　校正：.</b></p><p>　　3.三階龍格-庫(kù)塔方法【3】</p><p>　　類(lèi)似前面改進(jìn)的Euler方法公式的推導(dǎo)方法，將在處作Taylor展開(kāi)，然后再將在處作Taylor展開(kāi)，只要將兩個(gè)展開(kāi)式前四項(xiàng)相

9、同便有。于是得到三階龍格-庫(kù)塔公式為：</p><p>　　4.四階龍格-庫(kù)塔方法【3】 </p><p>　　類(lèi)似前面三階龍格-庫(kù)塔的推導(dǎo)方法，如果每步計(jì)算四次函數(shù)f（x，y）的值，完全類(lèi)似的，可以導(dǎo)出局部截?cái)嗾`差為的四階龍格-庫(kù)塔公式，其公式為：</p><p>　　四、程序代碼及運(yùn)算結(jié)果</p><p>　　用歐拉方法、改進(jìn)歐拉方法、3

10、階龍格－庫(kù)塔法以及4階龍格－庫(kù)塔法求解常微分方程初值問(wèn)題：</p><p>　　精確解為：y0=[ -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1050 -1.0000]</p><p><b>　　精確解圖形</b><

11、;/p><p><b>　　1、用歐拉法求解</b></p><p><b>　　程序如下：</b></p><p>　　建立函數(shù)文件cwfa1.m</p><p>　　function [x,y]=cwfa1(fun,x_span,y0,h)</p><p>　　x=x_spa

12、n(1):h:x_span(2);</p><p><b>　　y(1)=y0;</b></p><p>　　for n=1:length(x)-1 </p><p>　　y(n+1)=y(n)+h*feval(fun,x(n),y(n));</p><p><b>　　end</b></p&

13、gt;<p>　　x=x';y=y';</p><p>　　在MATLAB輸入以下程序：</p><p>　　>> clear all</p><p>　　>> fun=inline(' y^2*x ');</p><p>　　>> [x,y]=cwfa1(fu

14、n,[0,1],-2,0.1);</p><p><b>　　>> [x,y]</b></p><p>　　>> plot(x,y,'r+-')</p><p><b>　　ans =</b></p><p>　　0 -2.0000</p>

15、<p>　　0.1000 -2.0000</p><p>　　0.2000 -1.9600</p><p>　　0.3000 -1.8832</p><p>　　0.4000 -1.7768</p><p>　　0.5000 -1.6505</p><p>　　0.6000 -1.51

16、43</p><p>　　0.7000 -1.3767</p><p>　　0.8000 -1.2440</p><p>　　0.9000 -1.1202</p><p>　　1.0000 -1.0073</p><p><b>　　結(jié)果及其圖象：</b></p>&

17、lt;p>　　圖1 歐拉法解 </p><p>　　2、用改進(jìn)歐拉法求解：</p><p><b>　　程序如下：</b></p><p>　　建立函數(shù)文件cwfa2.m</p><p>　　function [x,y]=cwfa2(fun,x_span,y0,h)</p>

18、<p>　　x=x_span(1):h:x_span(2);</p><p><b>　　y(1)=y0;</b></p><p>　　for n=1:length(x)-1</p><p>　　k1=feval(fun,x(n),y(n));</p><p>　　y(n+1)=y(n)+h*k1;</

19、p><p>　　k2=feval(fun,x(n+1),y(n+1));</p><p>　　y(n+1)=y(n)+h*(k1+k2)/2;</p><p><b>　　end</b></p><p>　　x=x';y=y';</p><p>　　在MATLAB輸入以下程序：<

20、/p><p>　　>> clear all</p><p>　　>> fun=inline(' y^2*x'); </p><p>　　>> [x,y]=cwfa2(fun,[0,1],-2,0.1);</p><p><b>　　>> [x,y]</b><

21、;/p><p>　　>> plot(x,y,'r+-') </p><p><b>　　ans =</b></p><p>　　0 -2.0000</p><p>　　0.1000 -1.9800</p><p>　　0.2000 -1.9227</p&g

22、t;<p>　　0.3000 -1.8345</p><p>　　0.4000 -1.7239</p><p>　　0.5000 -1.6001</p><p>　　0.6000 -1.4711</p><p>　　0.7000 -1.3432</p><p>　　0.8000 -

23、1.2208</p><p>　　0.9000 -1.1066</p><p>　　1.0000 -1.0018</p><p><b>　　結(jié)果及其圖像：</b></p><p><b>　　圖3改進(jìn)歐拉解</b></p><p>　　3、用3階龍格—庫(kù)塔求解<

24、;/p><p><b>　　程序如下：</b></p><p>　　建立函數(shù)文件cwfa4.m</p><p>　　function [x,y]=cwfa4(fun,x_span,y0,h)</p><p>　　x=x_span(1):h:x_span(2);</p><p><b>　　y

25、(1)=y0;</b></p><p>　　for n=1:length(x)-1</p><p>　　k1=feval(fun,x(n),y(n));</p><p>　　k2=feval(fun,x(n)+h/2,y(n)+h/2*k1);</p><p>　　k3=feval(fun,x(n)+h,y(n)+h*(2*k2-

26、k1));</p><p>　　y(n+1)=y(n)+h*(k1+4*k2+k3)/6;</p><p><b>　　end</b></p><p>　　x=x';y=y';</p><p>　　在Matlab輸入以下程序：</p><p>　　clear all;</p&

27、gt;<p>　　fun=inline(' y^2*x');</p><p>　　[x,y]=cwfa4(fun,[0,1],-2 ,0.1);</p><p><b>　　[x,y] </b></p><p>　　plot(x,y, 'b*-')</p><p><b

28、>　　ans =</b></p><p>　　0 -2.0000</p><p>　　0.1000 -1.9803</p><p>　　0.2000 -1.9232</p><p>　　0.3000 -1.8350</p><p>　　0.4000 -1.7243</p>

29、;<p>　　0.5000 -1.6002</p><p>　　0.6000 -1.4707</p><p>　　0.7000 -1.3424</p><p>　　0.8000 -1.2196</p><p>　　0.9000 -1.1050</p><p>　　1.0000 -1

30、.0000</p><p><b>　　結(jié)果及圖像：</b></p><p>　　圖3 3階龍格-庫(kù)塔解</p><p>　　4、 用4階龍格—庫(kù)塔求解</p><p><b>　　程序如下：</b></p><p>　　建立函數(shù)文件cwfa3.m</p>&l

31、t;p>　　function [x,y]=cwfa3(fun,x_span,y0,h)</p><p>　　x=x_span(1):h:x_span(2);</p><p><b>　　y(1)=y0;</b></p><p>　　for n=1:length(x)-1</p><p>　　k1=feval(fun

32、,x(n),y(n));</p><p>　　k2=feval(fun,x(n)+h/2,y(n)+h/2*k1);</p><p>　　k3=feval(fun,x(n)+h/2,y(n)+h/2*k2);</p><p>　　k4=feval(fun,x(n+1),y(n)+h*k3);</p><p>　　y(n+1)=y(n)+h*(

33、k1+2*k2+2*k3+k4)/6;</p><p><b>　　end</b></p><p>　　x=x';y=y';</p><p>　　在MATLAB輸入以下程序：</p><p>　　>> clear all;</p><p>　　>> fun=

34、inline(' y^2*x ');</p><p>　　>> [x,y]=cwfa3(fun,[0,1],-2 ,0.1);</p><p><b>　　>> [x,y] </b></p><p>　　>> plot(x,y, 'b*-') </p><p

35、><b>　　ans =</b></p><p>　　0 -2.0000</p><p>　　0.1000 -1.9802</p><p>　　0.2000 -1.9231</p><p>　　0.3000 -1.8349</p><p>　　0.4000 -1.7241

36、</p><p>　　0.5000 -1.6000</p><p>　　0.6000 -1.4706</p><p>　　0.7000 -1.3423</p><p>　　0.8000 -1.2195</p><p>　　0.9000 -1.1050</p><p>　　1.

37、0000 -1.0000</p><p><b>　　結(jié)果及其圖象：</b></p><p>　　圖4 4階龍格-庫(kù)塔解</p><p>　　6、歐拉方法與改進(jìn)歐拉方法、3階龍格－庫(kù)塔法以及4階龍格－庫(kù)塔法得出的解得比較。</p><p>　　x=0:0.1:1;</p><p>　　y0=[

38、 -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1050 -1.0000]</p><p>　　y1=[ -2.0000 -2.0000 -1.9600 -1.8832 -1.7768 -1.6505 -1.5143 -1.3767 -1.

39、2440 -1.1202 -1.0073];</p><p>　　y2=[ -2.0000 -1.9800 -1.9227 -1.8345 -1.7239 -1.6001 -1.4711 -1.3432 -1.2208 -1.1066 -1.0018];</p><p>　　y3=[ -2.0000 -1.9803 -1.9232

40、 -1.8350 -1.7243 -1.6002 -1.4707 -1.3424 -1.2196 -1.1050 -1.0000];</p><p>　　y4=[ -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1005 -1.1000];<

41、;/p><p>　　plot(x,y0,'r*-')</p><p>　　hold on,plot(x,y1,'b-'),title('精確解與歐拉方法比較')</p><p>　　plot(x,y0,'r+-')</p><p>　　hold on,plot(x,y2,'b

42、-'),title('精確解與改進(jìn)歐拉方法比較')</p><p>　　plot(x,y0,'r+-')</p><p>　　hold on,plot(x,y3,'b-'),title('精確解與3階龍格—庫(kù)塔方法比較')</p><p>　　plot(x,y0,'r+-')&l

43、t;/p><p>　　hold on,plot(x,y4,'b-'),title('精確解與4階龍格—庫(kù)塔方法比較') </p><p>　　圖5 精確解與歐拉法比較</p><p>　　圖6精確解與改進(jìn)歐拉法比較</p><p>　　圖7精確解與3階-龍格庫(kù)塔方法比較</p><p>　

44、　圖8 精確解與4階龍格-庫(kù)塔方法比較</p><p>　　五、數(shù)值分析設(shè)計(jì)的gui界面</p><p><b>　　六、結(jié)果分析</b></p><p>　　由以上的各種方法與精確解的比較圖可以看出在計(jì)算精度上，四階龍格-庫(kù)塔方法的誤差最小，其次是三階龍格-庫(kù)塔方法，再次之就是改進(jìn)歐拉方法，歐拉方法誤差則比較大，所以四階龍格-庫(kù)塔方法得到最佳

45、的精度。而在計(jì)算量上面，相應(yīng)地，很明顯的四階龍格-庫(kù)塔方法也是最大，三階龍格-庫(kù)塔次之，其后為改進(jìn)歐拉方法，歐拉方法計(jì)算量最小。這樣的結(jié)果，說(shuō)明了運(yùn)用以上三種方法時(shí)，其計(jì)算量的多少與精度的大小成正比。我們?cè)趯?shí)際運(yùn)用與操作中，可以根據(jù)實(shí)際情況，選擇這4種方法中的其中一種最適合的，追求精度的話(huà)，可以使用三階或四階龍格-庫(kù)塔方法；而改進(jìn)的歐拉方法，在精度上和計(jì)算量上都表現(xiàn)得很出色，能夠滿(mǎn)足一般情況；而歐拉方法更主要的是適用于對(duì)的估計(jì)上，相應(yīng)的

46、，精度則有所欠缺。以上的選擇，都取決于具體的情況。</p><p><b>　　七、設(shè)計(jì)心得</b></p><p>　　這次實(shí)驗(yàn)花了較多的時(shí)間，先是選擇方程方面遇到了困難：有些方程在運(yùn)行過(guò)程中會(huì)不斷出錯(cuò)，錯(cuò)誤的原因有很多，有的是在輸入格式上有錯(cuò)，有的是在步長(zhǎng)和初值的取法上出錯(cuò)。在試了十多個(gè)方程之后終于找到一個(gè)自己滿(mǎn)意的。接下來(lái)的運(yùn)行中也出現(xiàn)了一些錯(cuò)誤。在用第一個(gè)方程

47、求精確解時(shí)由于沒(méi)有區(qū)分點(diǎn)乘和乘導(dǎo)致錯(cuò)誤，有些方程精確解時(shí)可以運(yùn)行但是在用歐拉方程和改進(jìn)的歐拉方法時(shí)就會(huì)報(bào)錯(cuò)，有的方程即使求的出來(lái)結(jié)果與精確解的誤差是很大的，顯然是錯(cuò)誤的。還有就是剛開(kāi)始犯了一個(gè)低級(jí)錯(cuò)誤：把原函數(shù)直接代入歐拉求解這樣當(dāng)然是不正確的。取初值對(duì)函數(shù)圖像的影響很大，如果初值取的不合適誤差也會(huì)非常大，所以要合理的選擇初值，在幾種求解方法中只有龍格—庫(kù)塔的求解與精確值最為接近幾乎是重合的。通過(guò)的后面的不同步長(zhǎng)求解的比較發(fā)現(xiàn)步長(zhǎng)越小與

48、精確值越接近?？傊谧鲞@次實(shí)驗(yàn)中對(duì)學(xué)習(xí)數(shù)學(xué)有了更多的體會(huì)，那就是一定要認(rèn)真，不能絲毫馬虎，一個(gè)小小的錯(cuò)誤就會(huì)導(dǎo)致所有的數(shù)據(jù)不能運(yùn)行。</p><p><b>　　八、參考文獻(xiàn)</b></p><p>　　[1] 劉衛(wèi)國(guó).《MATLAB程序設(shè)計(jì)與應(yīng)用程》. 北京: 高等教育出版社，2009</p><p>　　[2]奚如成 《數(shù)值分析方法》 中國(guó)