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

1、攵算H夬駕EAST CHINA INSTITUTE OF TECHNOLOGY數(shù)值分析課程設計題目：常微分方程的數(shù)值解法組員：鄧全飛 201820320185許曦 2018203201熊魯平20120320203余佳明 201820320206專業(yè)：信息與計算科學指導老師：劉唐偉、摘要3二、設計目的3三、理論基礎31. 歐拉公式【3】32. 改進Euler方法【3】33. 三階龍格-庫塔方法【3】34. 四階龍格-庫塔方法【3】4四、程序代碼及運算結果41、用歐拉法求解 52、 用改進歐拉法求解：63、用3階龍格一庫塔求解 74、用4階龍格一庫塔求解 9106、歐拉方法與改進歐拉方法、3階龍格

2、-庫塔法以及 4階龍格-庫塔法得岀的解得比較。五、數(shù)值分析設計的GUI界面13六、結果分析14七、設計心得14八、參考文獻14一、摘要在matlab環(huán)境下熟悉的運用計算機編程語言并結合龍格-庫塔法的理論 基礎對常微分方程初值問題進行求解，在運行完程序后以及對運行結果做出各 方面的分析和比較。二、設計目的用熟悉的計算機語言編程上機完成用歐拉方法、改進歐拉方法、3階龍格-庫塔法以及4階龍格-庫塔法求解常微分方程初值問題。三、理論基礎1.歐拉公式【3】作Taylor展開，得在點 將，那么當h充分小時，略去誤差項上J| ，用 近似替代廠1、廠近似替代廠1，并注意到，便得上述方法稱為Euler方法。2.

3、改進Euler方法【3】在應用梯形方法的迭代公式進行運算時，每迭代一次都要重新計算函數(shù)的值，且還要判斷何時可以終止或轉下一步計算。為了控制計算量和簡化算法，通常只迭代一次就轉入下一步計算。具體說，我們先用Euler公式求得一個初步的近似值，稱之為預測值，然后用梯形方法的迭代公式作一次迭代得 ，即將校正一次，這樣建立的預測-校正方法稱之為改進的Euler方法:預測：校正：.二-u fl3.三階龍格-庫塔方法【3】類似前面改進的Euler方法公式的推導方法，將在LI處作Taylor展開，然后再將 I在 I處作Taylor展開，只要將兩個展開式前四項相 同便有 一 =H。于是得到三階龍格-庫塔公式為

4、：【3】4.四階龍格-庫塔方法類似前面三階龍格-庫塔的推導方法，如果 每步計算四次函數(shù) f x=x_span(1:h:x_span(2 。 y(1=y0 。for n=1:length(x-1y(n+1=y(n+h*feval(fun,x(n,y(n end x=x 。 y=y 。在MATLA輸入以下程序： clear all fun=inlin e( yA2*x 。 x,y=cwfa1(fun,0,1,-2,0.1 x,y plot(x,y,r+-ans =0 -2.00000.1000-2.00000.2000-1.96000.3000-1.88320.4000-1.77680.5000-

5、1.65050.6000-1.51430.7000-1.37670.8000-1.24400.9000-1.12021.0000-1.0073結果及其圖象：圖1歐拉法解2、用改進歐拉法求解： 程序如下：建立函數(shù)文件cwfa2.mfun cti on x,y=cwfa2(fu n, x_spa n,y 0,h x=x_span( 1:h:x_span(2 。y(1=y0 ofor n=1:le ngth(x-1 k1=feval(f un,x(n,y(n。y(n+1=y(n+h*k1。k2=feval(fu n,x( n+1,y( n+1 y( n+1=y( n+h*(k1+k2/2。endx=

6、x。y=y。在MATLA輸入以下程序： clear all fun=inlin e( yA2*x。 x,y=cwfa2(fu n,0,1,-2,0.1 x,y plot(x,y,r+-ans =6 / 140.1000-1.98000.2000-1.92270.3000-1.83450.4000-1.72390.5000-1.60010.6000-1.47110.7000-1.34320.8000-1.22080.9000-1.10661.0000 -1.0018結果及其圖像：圖3改進歐拉解3、用3階龍格一庫塔求解程序如下：建立函數(shù)文件cwfa4.mfun cti on x,y=cwfa4(f

7、u n, x_spa n,y 0,h x=x_span( 1:h:x_span(2 。y(1=y0 ofor n=1:le ngth(x-1k1=feval(f un,x(n,y(n。k2=feval(fu n,x( n+h/2,y( n+h/2*k1 k3=feval(fu n,x (n +h,y( n +h*(2*k2-k1 y(n+1=y( n+h*(k1+4*k2+k3/6。endx=x。y=y。在Matlab輸入以下程序： clear all 。fun=i nli ne( yA2*x。x,y=cwfa4(fu n,0,1,-2 ,0.1 x,yplot(x,y, b*-ans =0

8、-2.00000.1000 -1.98030.2000 -1.92320.3000 -1.83500.4000 -1.72430.5000 -1.60020.6000 -1.47070.7000 -1.34240.8000 -1.21960.9000 -1.10501.0000 -1.0000結果及圖像:圖3 3階龍格-庫塔解4、用4階龍格一庫塔求解程序如下：建立函數(shù)文件 cwfa3.mfunction x,y=cwfa3(fun,x_span,y0,h x=x_span(1:h:x_span(2 。y(1=y0 。for n=1:length(x-1 k1=feval(fun,x(n,y(n

9、 。 k2=feval(fun,x(n+h/2,y(n+h/2*k1 k3=feval(fun,x(n+h/2,y(n+h/2*k2 k4=feval(fun,x(n+1,y(n+h*k3 y(n+1=y(n+h*(k1+2*k2+2*k3+k4/6 endx=x 。 y=y 。在MATLA輸入以下程序： clear all 。 fun=inlin e( yA2*x 。 x,y=cwfa3(fun,0,1,-2 ,0.1 x,y plot(x,y, b*-ans =0 -2.00000.1000-1.98020.2000-1.92310.3000-1.83490.4000-1.72410.50

10、00-1.60000.6000-1.47060.7000-1.34230.8000-1.21950.9000-1.10501.0000-1.0000結果及其圖象：圖4 4階龍格-庫塔解6、歐拉方法與改進歐拉方法、3階龍格-庫塔法以及 4階龍格-庫塔法得出的解得比較x=0:0.1:1 。yO=-2.0000-1.9802-1.9231-1.8349-1.7241-1.6000-1.4706 -1.3423 -1.2195 -1.1050 -1.0000y1= -2.0000-2.0000-1.9600-1.8832-1.7768-1.6505-1.5143 -1.3767 -1.2440 -1.

11、1202 -1.0073。y2= -2.0000-1.9800-1.9227-1.8345-1.7239-1.6001- 1.4711 -1.3432 -1.2208 -1.1066 -1.0018y3= -2.0000 -1.9803 -1.9232 -1.8350 -1.7243 -1.6002 - 1.4707 -1.3424 -1.2196 -1.1050 -1.0000。y4= -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1005 -1.1000。plot(x,y0,r*-精確解

12、與歐拉方法比較精確解與改進歐拉方法比較精確解與3階龍格一庫塔方法比較hold on ,plot(x,y1,b-,title( plot(x,y0,r+-hold on ,plot(x,y2,b-,title(plot(x,y0,r+-hold on ,plot(x,y3,b-,title(plot(x,y0,r+-精確解與4階龍格一庫塔方法比較hold on ,plot(x,y4,b-,title(圖5精確解與歐拉法比較圖6精確解與改進歐拉法比較圖7精確解與3階-龍格庫塔方法比較圖8精確解與4階龍格-庫塔方法比較五、數(shù)值分析設計的gui界面六、結果分析由以上的各種方法與精確解的比較圖可以看出在

13、計算精度上，四階龍格-庫塔方法的誤差最小，其次是三階龍格-庫塔方法，再次之就是改進歐拉方法，歐 拉方法誤差則比較大，所以四階龍格-庫塔方法得到最佳的精度。而在計算量上 面，相應地，很明顯的四階龍格-庫塔方法也是最大，三階龍格-庫塔次之，其 后為改進歐拉方法，歐拉方法計算量最小。這樣的結果，說明了運用以上三種 方法時，其計算量的多少與精度的大小成正比。我們在實際運用與操作中，可 以根據(jù)實際情況，選擇這 4種方法中的其中一種最適合的，追求精度的話，可 以使用三階或四階龍格-庫塔方法；而改進的歐拉方法，在精度上和計算量上都 表現(xiàn)得很出色，能夠滿足一般情況；而歐拉方法更主要的是適用于對的估計上，相應的，精度則有所欠缺。以上的選擇，都取決于具體的情況。七、設計心得這次實驗花了較多的時間，先是選擇方程方面遇到了困難：有些方程在運行過 程中會不斷出錯，錯誤的原因有很多，有的是在輸入格式上有錯，有的是在步 長和初值的取法上出錯。在試了十多個方程之后終于找到一個自己滿意的。接 下來的運行中也出現(xiàn)了一些錯誤。在用第一個方程求精確解時由于沒有區(qū)分點 乘和乘導致錯誤，有些方程精確解時可以運行但是在用歐拉方程和改進的歐拉 方法時就會報錯，有的方程即使求的出來結果與精確解的誤差是很大的，顯然 是錯誤的。還有就是剛開始犯了一個低級錯誤：把原函數(shù)直接代入歐拉求解這 樣當然是不正確的。取初值對函數(shù)圖像的影響很大，