课程设计报告-四阶Runge-Kutta方法

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方法的误差取决于步长的选取，因此，在实验的时候我们需要慎