數(shù)值計(jì)算課程設(shè)計(jì) 四階RungeKutta方法

1、 湖南工業(yè)大學(xué)課 程 設(shè) 計(jì)資 料 袋 理 學(xué)院（系、部） 2013 學(xué)年第 2 學(xué)期 課程名稱 數(shù)值計(jì)算方法 指導(dǎo)教師 職稱 副教授 學(xué)生姓名 專業(yè)班級(jí) 信息與計(jì)算科學(xué)班 學(xué)號(hào)學(xué)生姓名 專業(yè)班級(jí) 信息與計(jì)算科學(xué)1002班 學(xué)號(hào)學(xué)生姓名 專業(yè)班級(jí) 信息與計(jì)算科學(xué) 學(xué)號(hào) 題 目 四階Runge-Kutta方法 成 績(jī) 起止日期 2013 年 6 月24日 2013 年 7月 5日目 錄 清 單序號(hào)材 料 名 稱資料數(shù)量備 注1課程設(shè)計(jì)任務(wù)書(shū)12課程設(shè)計(jì)說(shuō)明書(shū)13張 湖南工業(yè)大學(xué)課程設(shè)計(jì)任務(wù)書(shū) 2012 2013 學(xué)年第 2 學(xué)期 理 學(xué)院（系、部） 信息與計(jì)算科學(xué) 專業(yè) 1002 班級(jí)課程名稱：

2、 數(shù)值計(jì)算方法 設(shè)計(jì)題目： 四階Runge-Kutta方法 完成期限：自 2013 年 6 月 24 日至 2013 年 7月 5 日共 2 周內(nèi)容及任務(wù)1、 設(shè)計(jì)題目：四階Runge-Kutta方法的應(yīng)用2、 設(shè)計(jì)目的：編寫關(guān)于四階Runge-Kutta Matlab程序求解微分方程的初值問(wèn)題。進(jìn)度安排起止日期工作內(nèi)容6.24 -6.26進(jìn)行選題及審題6.27 -6.30資料準(zhǔn)備并進(jìn)行計(jì)算 7.01 -7.05課程設(shè)計(jì)報(bào)告書(shū)些階段主要參考資料數(shù)值計(jì)算方法 黃云清 舒適編著 科學(xué)出版社指導(dǎo)教師（簽字）： 年 月 日系（教研室）主任（簽字）： 年 月 日數(shù)值計(jì)算方法設(shè)計(jì)說(shuō)明書(shū)四階Runge-Ku

3、tta方法起止日期： 2013 年 6 月 24 日 至 2013 年 7月 5 日學(xué)生姓名班級(jí)信息與計(jì)算科學(xué)班學(xué)號(hào) 成績(jī)指導(dǎo)教師(簽字)理學(xué)院（院、部）2013年7月5日目 錄一、 摘要5二、 問(wèn)題重述5三、 方法原理及實(shí)現(xiàn)5四、 計(jì)算公式或算法5五、 Matlab程序6六、 測(cè)試數(shù)據(jù)及結(jié)果6七、 結(jié)果分析10八、方法改進(jìn)10九、心得體會(huì)10十、參考文獻(xiàn)101、 摘要本課程設(shè)計(jì)主要內(nèi)容是用四階Runge-Kutta方法解決常微分方程組初值問(wèn)題的數(shù)值解法，通過(guò)分析給定題目使用Matlab編寫程序計(jì)算結(jié)果并繪圖，最后對(duì)計(jì)算結(jié)果進(jìn)行分析，得到結(jié)論。2、 問(wèn)題重述 在計(jì)算機(jī)上實(shí)現(xiàn)用四階Runge-K

4、utta求一階常微分方程初值問(wèn)題 的數(shù)值解，并利用最后繪制的圖形直觀分析近似解與準(zhǔn)確解之間的比較。三、方法原理及實(shí)現(xiàn)龍格-庫(kù)塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對(duì)誤差進(jìn)行抑制，所以其實(shí)現(xiàn)原理也較復(fù)雜。該算法是構(gòu)建在數(shù)學(xué)支持的基礎(chǔ)之上的。龍格庫(kù)塔方法的理論基礎(chǔ)來(lái)源于泰勒公式和使用斜率近似表達(dá)微分，它在積分區(qū)間多預(yù)計(jì)算出幾個(gè)點(diǎn)的斜率，然后進(jìn)行加權(quán)平均，用做下一點(diǎn)的依據(jù)，從而構(gòu)造出了精度更高的數(shù)值積分計(jì)算方法。如果預(yù)先求兩個(gè)點(diǎn)的斜率就是二階龍格庫(kù)塔法，如果預(yù)先取四個(gè)點(diǎn)就是四階龍格庫(kù)塔法。經(jīng)典的方法是一個(gè)四階的方法，它的計(jì)算公式是：四、計(jì)算

5、公式或算法1 輸入（編寫或調(diào)用計(jì)算的函數(shù)文件），2 3For End 4輸出五、Matlab 程序x=a:h:b;y(1)=y1;n=(b-a)/h+1;for i=2:n fk1=f(x(i-1),y(i-1); fk2=f(x(i-1)+h/2,y(i-1)+fk1*h/2); fk3=f(x(i-1)+h/2,y(i-1)+fk2*h/2); fk4=f(x(i-1)+h,y(i-1)+fk3*h); y(i)=y(i-1)+h*(fk1+2*fk2+2*fk3+fk4)/6;endy六、測(cè)試數(shù)據(jù)及結(jié)果用調(diào)試好的程序解決如下問(wèn)題：應(yīng)用經(jīng)典的四階Runge-Kutta方法解初值問(wèn)題 ?。?

6、） 步驟一：編寫函數(shù)具體程序.1.求解解析解程序：dsolve(Dy=(y2+y)/t,y(1)=-2,t)結(jié)果：2.綜合編寫程序如下：a=1;b=3;h=0.5;y(1)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;for i=2:n k1=(y(i-1)2+y(i-1)/x(i-1); k2=(y(i-1)+h*k1/2)2+(y(i-1)+h*k1/2)/(x(i-1)+h/2); k3=(y(i-1)+h*k2/2)2+(y(i-1)+h*k2/2)/(x(i-1)+h/2); k4=(y(i-1)+h*k3)2+(y(i-1)+h*k3)/(x(i-1)+h); y

7、(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;% 四階Runge-Kutta公式解 x(i)=x(i-1)+h; %有解區(qū)間的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %誤差項(xiàng)endx y yy s（2）步驟二：執(zhí)行上述Runge-Kutta算法，計(jì)算結(jié)果為1.00001.50002.00002.50003.0000-2.0000-1.4954-1.3306-1.2480-1.1985-2.0000-1.5000-1.3333-1.2500-1.200000.00460.00280.00200.0015（3）使用

8、Matlab繪圖函數(shù)“plot(x,y)”繪制問(wèn)題數(shù)值解和解析解的圖形。數(shù)值解的圖形：plot(x,y)解析解的圖形plot(x,yy)（4）使用Matlab中的ode45求解，并繪圖。編寫函數(shù)如下：%ode.mfunction dy=ode(x,y)dy=(y2+y)/x; T,Y=ode45(ode,1 3,-2);plot(T,Y)運(yùn)行結(jié)果如下：7、 結(jié)果分析由圖可知此方法與精確解的契合度非常好，基本上與精度解保持一致，由此可見(jiàn)四階Runge-Kutta方法是一種高精度的單步方法。8、 方法改進(jìn)同時(shí)，由于誤差的存在，我們總想盡可能的是誤差趨近于零，常用的就是傳統(tǒng)的增加取值的個(gè)數(shù)。最后，我

9、們通過(guò)改變步長(zhǎng)來(lái)進(jìn)行改進(jìn)。具體實(shí)現(xiàn)：（1）h=0.1a=1;b=3;h=0.1;y(1)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;for i=2:n k1=(y(i-1)2+y(i-1)/x(i-1); k2=(y(i-1)+h*k1/2)2+(y(i-1)+h*k1/2)/(x(i-1)+h/2); k3=(y(i-1)+h*k2/2)2+(y(i-1)+h*k2/2)/(x(i-1)+h/2); k4=(y(i-1)+h*k3)2+(y(i-1)+h*k3)/(x(i-1)+h); y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;% 四階Runge

10、-Kutta公式解 x(i)=x(i-1)+h; %有解區(qū)間的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %誤差項(xiàng)endx y yy s結(jié)果：（2） h=0.2a=1;b=3;h=0.2;y(1)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;for i=2:n k1=(y(i-1)2+y(i-1)/x(i-1); k2=(y(i-1)+h*k1/2)2+(y(i-1)+h*k1/2)/(x(i-1)+h/2); k3=(y(i-1)+h*k2/2)2+(y(i-1)+h*k2/2)/(x(i-1)+h/2); k4=

11、(y(i-1)+h*k3)2+(y(i-1)+h*k3)/(x(i-1)+h); y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;% 四階Runge-Kutta公式解 x(i)=x(i-1)+h; %有解區(qū)間的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %誤差項(xiàng)endx y yy s結(jié)果： （3） h=0.4a=1;b=3;h=0.4;y(1)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;for i=2:n k1=(y(i-1)2+y(i-1)/x(i-1); k2=(y(i-1)+h*k1/2)

12、2+(y(i-1)+h*k1/2)/(x(i-1)+h/2); k3=(y(i-1)+h*k2/2)2+(y(i-1)+h*k2/2)/(x(i-1)+h/2); k4=(y(i-1)+h*k3)2+(y(i-1)+h*k3)/(x(i-1)+h); y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;% 四階Runge-Kutta公式解 x(i)=x(i-1)+h; %有解區(qū)間的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %誤差項(xiàng)endx y yy s結(jié)果：通過(guò)上述的一些結(jié)果得出，四階的Runge-Kutta方法的誤差取決于步長(zhǎng)的選取，因此，在實(shí)驗(yàn)的時(shí)候我們需要慎重的選取。一方面：我們要減少誤差，另一方面：我們也需要盡可能的減少計(jì)算次數(shù)。9、 心得體會(huì) 課程設(shè)計(jì)，至今我們小組三人感慨頗多，的確，從我們參考，設(shè)計(jì)到定稿，從理論到實(shí)踐，在整整兩星期的時(shí)間里，可以說(shuō)是苦多于甜，但是可以學(xué)到很多很多的東西，同時(shí)不僅可以鞏固以前所學(xué)過(guò)的知識(shí)，而且學(xué)到很多在書(shū)本上沒(méi)有學(xué)到的知識(shí)。通過(guò)