数值计算课程设计 四阶RungeKutta方法

1、 湖南工业大学课 程 设 计资 料 袋 理 学院（系、部） 2013 学年第 2 学期 课程名称 数值计算方法 指导教师 职称 副教授 学生姓名 专业班级 信息与计算科学班 学号学生姓名 专业班级 信息与计算科学1002班 学号学生姓名 专业班级 信息与计算科学 学号 题 目 四阶Runge-Kutta方法 成 绩 起止日期 2013 年 6 月24日 2013 年 7月 5日目 录 清 单序号材 料 名 称资料数量备 注1课程设计任务书12课程设计说明书13张 湖南工业大学课程设计任务书 2012 2013 学年第 2 学期 理 学院（系、部） 信息与计算科学 专业 1002 班级课程名称：

2、 数值计算方法 设计题目： 四阶Runge-Kutta方法 完成期限：自 2013 年 6 月 24 日至 2013 年 7月 5 日共 2 周内容及任务1、 设计题目：四阶Runge-Kutta方法的应用2、 设计目的：编写关于四阶Runge-Kutta Matlab程序求解微分方程的初值问题。进度安排起止日期工作内容6.24 -6.26进行选题及审题6.27 -6.30资料准备并进行计算 7.01 -7.05课程设计报告书些阶段主要参考资料数值计算方法 黄云清 舒适编著 科学出版社指导教师（签字）： 年 月 日系（教研室）主任（签字）： 年 月 日数值计算方法设计说明书四阶Runge-Ku

3、tta方法起止日期： 2013 年 6 月 24 日 至 2013 年 7月 5 日学生姓名班级信息与计算科学班学号 成绩指导教师(签字)理学院（院、部）2013年7月5日目 录一、 摘要5二、 问题重述5三、 方法原理及实现5四、 计算公式或算法5五、 Matlab程序6六、 测试数据及结果6七、 结果分析10八、方法改进10九、心得体会10十、参考文献101、 摘要本课程设计主要内容是用四阶Runge-Kutta方法解决常微分方程组初值问题的数值解法，通过分析给定题目使用Matlab编写程序计算结果并绘图，最后对计算结果进行分析，得到结论。2、 问题重述 在计算机上实现用四阶Runge-K

4、utta求一阶常微分方程初值问题 的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的比较。三、方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。经典的方法是一个四阶的方法，它的计算公式是：四、计算

5、公式或算法1 输入（编写或调用计算的函数文件），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六、测试数据及结果用调试好的程序解决如下问题：应用经典的四阶Runge-Kutta方法解初值问题 取（1

6、） 步骤一：编写函数具体程序.1.求解解析解程序：dsolve(Dy=(y2+y)/t,y(1)=-2,t)结果：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; %有解区间的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %误差项endx y yy s（2）步骤二：执行上述Runge-Kutta算法，计算结果为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绘图函数“plot(x,y)”绘制问题数值解和解析解的图形。数值解的图形：plot(x,y)解析解的图形plot(x,yy)（4）使用Matlab中的ode45求解，并绘图。编写函数如下：%ode.mfunction dy=ode(x,y)dy=(y2+y)/x; T,Y=ode45(ode,1 3,-2);plot(T,Y)运行结果如下：7、 结果分析由图可知此方法与精确解的契合度非常好，基本上与精度解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。8、 方法改进同时，由于误差的存在，我们总想尽可能的是误差趋近于零，常用的就是传统的增加取值的个数。最后，我

9、们通过改变步长来进行改进。具体实现：（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; %有解区间的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %误差项endx y yy s结果：（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; %有解区间的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %误差项endx y yy s结果： （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; %有解区间的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i); %误差项endx y yy s结果：通过上述的一些结果得出，四阶的Runge-Kutta方法的误差取决于步长的选取，因此，在实验的时候我们需要慎重的选取。一方面：我们要减少误差，另一方面：我们也需要尽可能的减少计算次数。9、 心得体会 课程设计，至今我们小组三人感慨颇多，的确，从我们参考，设计到定稿，从理论到实践，在整整两星期的时间里，可以说是苦多于甜，但是可以学到很多很多的东西，同时不仅可以巩固以前所学过的知识，而且学到很多在书本上没有学到的知识。通过