常微分方程數(shù)值解法課程設計報告

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

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

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

4、lt;/p><p><b>　　六、結果分析14</b></p><p><b>　　七、設計心得14</b></p><p><b>　　八、參考文獻14</b></p><p><b>　　一、摘要</b></p><p>　

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

6、lt;p><b>　　三、理論基礎</b></p><p><b>　　1.歐拉公式【3】</b></p><p>　　在點 將作Taylor展開，得</p><p>　　，那么當h充分小時，略去誤差項，用近似替代、近似替代，并注意到，便得</p><p>　　上述方法稱為Euler方法。&l

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

8、er方法：</p><p><b>　　預測：，</b></p><p><b>　　校正：.</b></p><p>　　3.三階龍格-庫塔方法【3】</p><p>　　類似前面改進的Euler方法公式的推導方法，將在處作Taylor展開，然后再將在處作Taylor展開，只要將兩個展開式前四項相

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

10、階龍格－庫塔法以及4階龍格－庫塔法求解常微分方程初值問題：</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>　　結果及其圖象：</b></p>&

17、lt;p>　　圖1 歐拉法解 </p><p>　　2、用改進歐拉法求解：</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>　　結果及其圖像：</b></p><p><b>　　圖3改進歐拉解</b></p><p>　　3、用3階龍格—庫塔求解<

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>　　結果及圖像：</b></p><p>　　圖3 3階龍格-庫塔解</p><p>　　4、 用4階龍格—庫塔求解</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>　　結果及其圖象：</b></p><p>　　圖4 4階龍格-庫塔解</p><p>　　6、歐拉方法與改進歐拉方法、3階龍格－庫塔法以及4階龍格－庫塔法得出的解得比較。</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('精確解與改進歐拉方法比較')</p><p>　　plot(x,y0,'r+-')</p><p>　　hold on,plot(x,y3,'b-'),title('精確解與3階龍格—庫塔方法比較')</p><p>　　plot(x,y0,'r+-')&l

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

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

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

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

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

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