课程设计报告-四阶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求一阶常微分方程初值问题；y'(x)

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

3、是：yn+=yn+h(Ki+2K2十2(十七)6Ki=f(xn,yn)一hh*K2=fM+,yn+3Ki)K3=f(Xnh,2Kh-2十ynK4-f(xnh,ynhKa)四、计算公式或算法Vo=yX01.输入a,b,n,y0,f(x,y)(编写或调用计算f(x,y)的函数文件F(x,y»,X0二a,Xn二bxi1-xoi-1hhKi=fXi工手K2=fXihi,xhiKiK3=fhi,yhiK2K4=fXyh,xhK3hyi=yiKi2K22K3K46End4.输出y1,y2yn.五、Matlab程序x=a:h:b;y(i)=yi;n=(b-a)/h+i;fori=2:nfki=f(

4、x(i-i),y(i-i);fk2=f(x(i-i)+h/2,y(i-i)+fki*h/2);fk3=f(x(i-i)+h/2,y(i-i)+fk2*h/2);fk4=f(x(i-i)+h,y(i-i)+fk3*h);y(i)=y(i-i)+h*(fki+2*fk2+2*fk3+fk4)/6;endy六、测试数据及结果用调试好的程序解决如下问题：应用经典的四阶Runge-Kutta方法解初值问题iny=-(y2y)1MtM3yt(yy),3,取h=0.5y(i)=-2,(i)步骤一：编写函数具体程序1 .求解解析解程序：dsolve('Dy=(yA2+y)/t','y(

5、1)=-2','t')结果：ans=2 .综合编写程序如下：a=1;b=3;h=0.5;y(i)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;fori=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(i-1)+h*(k1+2*k2+2*k3+k

6、4)/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算法，计算结果为X1.00001.50002.00002.50003.0000V-2.0000-1.4954-1.3306-1.2480-1.1985y(x)-2.0000-1.5000-1.3333-1.2500-1.2000s00.00460.00280.00200.00151.0000一2.0000-2.

7、000001.5000-L4954一1500G00462.odoo-L330S-L3333口.00282,5000-1.2480-1,25000.00203.0000-L1985-1.2000D.0015(3)使用Matlab绘图函数“plot(x,y)”绘制问题数值解和解析解的图形。数值解的图形:Plot(x,y)解析解的图形Plot(x,yy)(4)使用Matlab中的ode45求解，并绘图编写函数如下：%ode.mfunctiondy=ode(x,y)dy=(yA2+y)/x;T,Y=ode45('ode',13,-2);plot(T,Y)运行结果如下：七、结果分析由图可

8、知此方法与精确解的契合度非常好，基本上与精度解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。八、方法改进同时，由于误差的存在，我们总想尽可能的是误差趋近于零，常用的就是传统的增加取值的个数。最后，我们通过改变步长来进行改进。具体实现：(1)h=0.1a=1;b=3;h=0.1;y(i)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;fori=2:nk1=(y(i-1)A2+y(i-1)/x(i-1);k2=(y(i-1)+h*k1/242+(y(i-1)+h*k1/2)/(x(i-1)+h/2);k3=(y(i-1)+h*k2/2)A2+(y(i-1)

9、+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+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'结果：1.0000-2,0000-2.0000Q1.1000-LS333-1.83330.00001.2000-L71437.71430.00001.3000-k62

10、50-1.62500.00001.4D00-L5556-1.55560.00001.5000-1.5000-1,50000.00001.6000-1.4545-L45450.00001.7000-L4167-1.41670,00001.8000-1.3846-L38460,00001.9000-1.3571-L3E710.00002.0000-1.3333-1.33330.00002.10OQ-1.3125-1.31250.00002.2000-1.2941-1.29410.00002.3000-1.2778-1.27780,00002.4000-1.2632-1.26320.00002.50

11、00-1.2500-1.25000.00002.6000-1.2381-1.23810.00002.7100-1.2273-1.22730.00002.8000-L2174-1.21740.00002.9000-L2083-L20830.00003.0000-L2000-1.20000.0000h=0.2a=1;b=3;h=0.2;y(i)=-2;x(1)=a;n=(b-a)/h+1;yy(i)=-2;fori=2:nk1=(y(i-1)A2+y(i-1)/x(i-1);k2=(y(i-1)+h*k1/242+(y(i-1)+h*k1/2)/(x(i-1)+h/2);k3=(y(i-1)+h*

12、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(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=0100OQ.oooo0.00000.00000.00000.00000.00000.00000.00000,00001.

13、0000-2.0000-2.00001.20007.7142-1.71431.4000-L5E55-1.55561.6OCO-1.4545-1.45451.8000-1.3846-1.38462.0000-1.3333-1.33332,2000-1,2941-1.29412,4000-1.2631-1.26322.6000-L2381-1.23812.BODO-L21M-1.21743.0000-1,2000-1.2000(3)h=0.4a=1;b=3;h=0.4;y(1)=-2;x(1)=a;n=(b-a)/h+1;yy(1)=-2;fori=2:nk1=(y(i-1)A2+y(i-1)/x

14、(i-1);k2=(y(i-1)+h*k1/242+(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);四阶Runge-Kutta公式解y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;%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-L5640-1.55560.00161.8000-L3836-1.38460.00102,2000-1.2934-1.29410.0007Z6000-1.2375-1.23810.00063.0000-L1995-1.20000.0005通过上述的一些结果得出，四阶的Runge-Kutta方法的误差取决于步长的选取，因此，在实验的时候我们需要慎重的选取。一方面：我们要减少误差，另一方面：我们也需要尽可能的减少计算次数。九、心得体会