《数值计算课程设计-- 四阶Runge-Kutta方法》由会员分享，可在线阅读，更多相关《数值计算课程设计-- 四阶Runge-Kutta方法（13页珍藏版）》请在金锄头文库上搜索。
1、1湖南工业大学课 程 设 计资 料 袋理 学院（系、部） 2013 学年第 2 学期 课程名称 数值计算方法 指导教师 职称 副教授 学生姓名 专业班级 信息与计算科学班 学号学生姓名 专业班级 信息与计算科学 1002 班 学号学生姓名 专业班级 信息与计算科学 学号 题 目 四阶 Runge-Kutta 方法 成 绩 起止日期 2013 年 6 月 24 日 2013 年 7 月 5 日目 录 清 单序号材 料 名 称资料数量备 注1课程设计任务书12课程设计说明书13张湖南工业大学2课程设计任务书2012 2013 学年第学年第 2 学期学期理 学院（系、部） 信息与计算科学 专业 10
2、02 班级课程名称： 数值计算方法 设计题目： 四阶 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课程设计报告书些阶段 主 要 参 考 资 料数值计算方法 黄云清 舒适编著 科学出版社指导教师 （签字）： 年 月 日系（教研室）主任（
3、签字）： 年 月 日3数值计算方法设计说明书四阶 Runge-Kutta 方法起止日期： 2013 年 6 月 24 日 至 2013 年 7 月 5 日学生姓名班级信息与计算科学班学号成绩指导教师 (签字 )理学院（院、部）理学院（院、部）2013 年年 7 月月 5 日日4目目 录录一、 摘要5二、 问题重述5三、 方法原理及实现5四、 计算公式或算法5五、 Matlab 程序6六、 测试数据及结果6七、 结果分析510八、方法改进10九、心得体会10十、参考文献101 1、摘要摘要本课程设计主要内容是用四阶 Runge-Kutta 方法解决常微分方程组初值问题的数值解法，通过分析 给定题
4、目使用 Matlab 编写程序计算结果并绘图，最后对计算结果进行分析，得到结论。62 2、问题重述问题重述在计算机上实现用四阶 Runge-Kutta 求一阶常微分方程初值问题 1,yaybaxyxfxy的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的比较。三、方法原理及实现三、方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采 取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔 方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然 后进行加权
5、平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点 的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。 经典的方法是一个四阶的方法，它的计算公式是：RK112341213243(22)6 (,)(,)22(,)22 (,)nnnnnnnnnnhyyKKKKKf xyhhKf xyKhhKf xyKKf xhyhK 四、计算公式或算法四、计算公式或算法1 输入（编写或调用计算的函数文件） ，yxfynba,0yxf,yxF,2b-ahn（） 00xyy bxaxn ,03For ni:1 101111211111311112411312 ,iiiiiii
6、iixxihhhKf xyKf xh yh KKf xh yh KKf xh yhK743211226KKKKhyyiiEnd 4输出.,21nyyy 五、五、MatlabMatlab 程序程序x=a:h:b; y(1)=y1; n=(b-a)/h+1; for i=2:nfk1=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; end y六、测试
7、数据及结果六、测试数据及结果用调试好的程序解决如下问题： 应用经典的四阶 Runge-Kutta 方法解初值问题取21(),13,(1)2,yyytt y 0.5h=（1 1） 步骤一：编写函数具体程序步骤一：编写函数具体程序. . 1.1.求解解析解程序：求解解析解程序：dsolve(Dy=(y2+y)/t,y(1)=-2,t)结果：结果：2.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:n8k1=(y(i-1)2+y(i-1)/x(i-1);k2=(y(i-1)+h*
8、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-Kutta 公式解x(i)=x(i-1)+h; %有解区间的值yy(i)=-x(i)/(x(i)-1/2); %解析解s(i)=abs(y(i)-yy(i); %误差项end x y yy s（2 2）步骤二：执行上述）步骤二：执行上述 Runge-Ku
9、ttaRunge-Kutta 算法，计算结果为算法，计算结果为ix1.00001.50002.00002.50003.0000iy-2.0000-1.4954-1.3306-1.2480-1.1985 iy x-2.0000-1.5000-1.3333-1.2500-1.2000is00.00460.00280.00200.0015（3 3）使用）使用 MatlabMatlab 绘图函数绘图函数“plot(x,y)”“plot(x,y)”绘制问题数值解和解析解的图形。绘制问题数值解和解析解的图形。 数值解的图形：plot(x,y)9解析解的图形plot(x,yy)（4 4）使用）使用 Matl
10、abMatlab 中的中的 ode45ode45 求解，并绘图。求解，并绘图。10编写函数如下： %ode.mfunction dy=ode(x,y)dy=(y2+y)/x; T,Y=ode45(ode,1 3,-2);plot(T,Y)运行结果如下：7、结果分析结果分析由图可知此方法与精确解的契合度非常好，基本上与精度解保持一致，由此可见四阶 Runge-Kutta 方法是一种高精度的单步方法。8、方法改进方法改进同时，由于误差的存在，我们总想尽可能的是误差趋近于零，常用的就是传统的增加取值的个数。 最后，我们通过改变步长来进行改进。 具体实现： （1 1）h=0.1h=0.1a=1; b=
11、3; h=0.1; y(1)=-2; x(1)=a;11n=(b-a)/h+1; yy(1)=-2; for i=2:nk1=(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-Kutta 公式解x(i)=x(i-1)+h; %有解区
12、间的值yy(i)=-x(i)/(x(i)-1/2); %解析解s(i)=abs(y(i)-yy(i); %误差项end x y yy s结果：结果：（2 2）h=0.2h=0.2 a=1; b=3; h=0.2; y(1)=-2; x(1)=a;12n=(b-a)/h+1; yy(1)=-2; for i=2:nk1=(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-Kutta 公式解x(i)=x(i-1)+h; %有解区间的值yy(i)=-x(i)/(x(i)-1
