版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、计算机数值方法课程设计报告题目四阶Runge-Kutta 方法学生姓名班级计科12学号成绩指导教师延安大学计算机学院2014年9月1日目录一、摘要5二、问题重述5三、方法原理及实现5四、计算公式或算法5五、Matlab 程序6六、测试数据及结果6七、结果分析10八、方法改进10九、心得体会10十、参考文献10一、摘要本课程设计主要内容是用四阶Runge-Kutta方法解决常微分方程组初值问题的数值解法,首先分析题目内容和要求,然后使用Matlab编写程序计算结果并绘图,最后对计算结果进行分析并得出结论。二、问题描述在计算机上实现用四阶 Runge Kutta求一阶常微分方程初值问题Iy x f
2、 x, y ,x a,b,y a yi的数值解,并利用最后绘制的图形直观分析近似解与准确解之间的比较三、方法原理及实现龙格一库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法.由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂.该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法.如果预先求两个点的斜率就是二阶龙格库塔法,经典的R K方法是一个四阶的方法,它的计算公式是yn 1h/izyn(K162K2 2
3、K3K4)K1f(Xn, Yn)号KJK2hf (Xn 刁,ynK3hf (xn , ynK4f (xn h, ynhK3)如果预先取四个点就是四阶龙格库塔法。四、计算公式或算法1.输入 a,b,n, y0, fx, y (编写或调用计算f x, y的函数文件F x, y ),yoy xoxo a,Xn b3 . For i 1 : nx 1Xoi1 hh1h2K1fXi1, yi 1K2fXi1h1, yi1hKK3fXi1h1, yi1hK 2K4fXi1h, yi1 hK3hyiVi 1K12K2 2K3 K6End4 输出 yi, y2yn.五、Matlab程序x=a : h:b ;y
4、(i) =yi;n=(b-a)/h+1;for i=2:nfk仁f(x(i-1),y( i D);fk2=f(X(i-1)+h/2,y(i 1)+fk1 * h/2 )fk3=f(X(i1)+h/2,y (i-1)+fk2*h/2)fk4=f(X(i 1)+h,y(i1)+fk3*h );y(i ):=y( i-1)+h * (fk1+2*fk2+2*fk3+fk4)/6endy六、测试数据及结果用调试好的程序解决如下问题:应用经典的四阶 Runge Kutta方法解初值问题1 2y 1(yy), 1 t 3,取 h=0.5y(1)2,(1) 步骤一:编写函数具体程序1。求解解析解程序:dso
5、lve(Dy=(yA2+y)/t结果:ans =2.综合编写程序如下:a=1;b=3;h=0.5;y (1) =-2;x (1) =a;n= (ba) /h+1;yy(1 ) =2;for i=2: nk1=(y(i-1 ) A2+y(i 1)/x(i-1);k2=(y(i-1)+h*k1/2)A2+(y(i-1) +h* k1/2)/(x(i 1)+h/2 );k3=(y(i 1) +h* k2/2)A2+(y(i 1)+h*k2/2) )/ (x(i-1 )+h/2);k4=(y (i 1) +h* k3)A2+ (y (i-1)+h* k3) ) /(x(i-1) +h);y (i)=y
6、 (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算法,计算结果为Xi1.00001.50002.00002.50003.0000yi2.00001.49541.33061。 2480-1。1985y Xi-2 。 0000-1.5000-1。3333-1。25001。 2000s00.00460。
7、 00280.00200.00151.0000-2. 0000-2.000001.5000-1,9&4-L 5000- 00452.0000-L33330.00282.5COO-1_ 2480-1,2500D. 00203*ocoo1. 198S-L 20000. 0015(3) 使用Matlab绘图函数“ plot(x,y)绘制问题数值解和解析解的图形数值解的图形:plot(x , y)解析解的图形plot (x, yy)PxELire I(4) 使用Matlab中的ode45求解,并绘图.编写函数如下:%ode.mfunction dy=ode( x, y)dy=(yA2+y ) /x;T
8、, Y=ode45(ode , 1 3 , T);plot(T , Y)运行结果如下:七、结果分析由图可知此方法与精确解的契合度非常好, 基本上与精度解保持一致, 由此可见四阶 Runge-Kutta 方 法是一种高精度的单步方法。八、方法改进同时,由于误差的存在, 我们总想尽可能的是误差趋近于零 , 常用的就是传统的增加取值的个数。 最后, 我们通过改变步长来进行改进。具体实现:(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:nk1=(y(i-1) A2+y (i-1 ) )/x(i-1);k
9、2=(y (i-1)+h * k1/2F2+(y(i1)+h*k1/2 ) )/ (x(i 1)+h/2 );k3=(y(i 1 ) +hk2/2)A2+ (y(i-1 ) +h* k2/2) ) / ( x(i 1)+h/2 )k4=(y(i 1)+h*k3)A2+ (y(i-1)+h* k3) ) / (x(i-1)+h);y(i)=y(i-1) +h( k1+2k2+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
10、); % 误差项endx y yy s结果:ansI. 00002. 00002.000001.1000-1.8333-1.83330. aoooI.2000-1.1431.71430. ouoo1.3000-1_ 6250-L62500, oooc1.4000-1.55551.55560. 00001. 5000-1.5000-1-50000. ooaoI.eOQO-1.45451.45450. 00001. 7000-L4L6F-1-4LS70. oooo1.8000-L 3646-1.3846CL 00001. 9000-1.3571-1.351o. ooao2.0000-1.3333-
11、1-3333n oooo2.1000-1.3125-1.31250, oooo2.2000-1.2941-1-2941a. oooo2.3000-L 2778-1-2778o. ooao2,4000-L 2632-1-2032o. aooo2.5000-1.2B00-1.2500o. ooao2.6000-L2381-1.23810. 00002.FOOD-L 22F3-1.2273o.ooao2_0OOO-1_2174-1.2174oooo2.9000-L2083-1.2083o. oooo3.0000-1.2000-1.2000o. oooo(2)h=0。 2a=1;b=3;h=0.2;y
12、(1)=-2;x (1) =a;n= (b-a)/h+1 ;yy(1 ) =-2;for i=2: nk1 =(y (i 1)A2+y(i-1)/x(i -1);k2=(y (i -1)+h * k1/2 )A2+ (y(i 1)+h * k1/2)/(x(i-1 ) +h/2);k3=(y (i-1)+h*k2/2 )A2+ (y(i-1)+h*k2/2)/(x(i 1)+h/2 );k4=(y (i-1)+h* k3) A2+ (y(i1)+h * k3) )/(x(i-1)+h);y (i ) =y(i 1)+h * (k1+2 * k2+2 * k3+k4) /6 ;% 四阶 Rung
13、e-Kutta 公式解yy(i)= x(i)/(x(i )-1/2 );%解析解s(i ) =abs (y (i)-yy(i);%误差项endxy yy s 结果:arts =1. 0000-2.0D00-2.000001. 2000-1.7142-1. 71430.0000L4000-1.5555-1.55500.00001-6000-1.4545-L45450.00001.800Q-1.3845-1.3846o.aaoo2. 0000-1.3333-1-33330. 00002. 2000-1.2941-1.2941o.oooa2.4000-1.2631-1.26320.00002. 00
14、00-1.2381-1.23810.0000N 8000-1.2174-1.21740.00003-0000-1.2000-L20D00.0000(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(i-1 );(y (i 1) +h* k1/2 ) / (x(i 1)+h/2);k仁(y(i 1)A2+y (i-1)/xk2=(y(i-1 ) +h* k1/2F2+k3=(y(i 1) +h*k2/2 ) A2+(y(i-l)+h*k2/2) )/(x(i 1) +h/2 )k4=(y(i-1
15、) +h* k3) A2+(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结果:ans =1. 0000-2. 0000-2. 000001.4000-1.E640-1.EE560.ooie1.8000-1.3836-L 38400. 00102.2000-L2934-1.29410. 00072,6000-L23F5-1.23810. 00063.QQOO-L 1996-1.20000.0005通过上述的一些结果得出,四阶的Runge Kutta方法的误差取决于步长的选取,因此,在实验的时候我们需要慎
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 2024届广东肇庆市高三4月质量调研(二模)考试数学试题
- 餐饮店合同补充协议范本
- 财产处份协议书
- 亳州公证处合同公证收费标准
- 北京市租房标准合同
- 山西省2024八年级物理上册第三章物态变化第2节熔化和凝固第2课时熔化和凝固的应用课件新版新人教版
- 设备维修班长述职报告
- 湖南省益阳市赫山区箴言龙光桥学校2024-2025学年四年级上学期期中考试数学试题(无答案)
- 《J类船用筒形观察器》
- 广西柳州市2024-2025学年七年级上学期11月期中考试数学试题(含答案)
- 公路超限超载检测站设计指南(试点工程版)
- 化工生产开停车方案
- 高精度时间同步及定位技术应用白皮书
- 小学科学教育科学三年级上册空气《风的成因》教案
- DB13 2863-2018 炼焦化学工业大气污染物超低排放标准
- 火炬系统水封罐计算
- 怎样写好文学短评课件(15张PPT)
- 医院卒中中心建设方案
- 大学新生社团招新宣传PPT模板课件
- T∕CCMSA 10411-2020 铸铝门
- 氧环境下材料的选择
评论
0/150
提交评论