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

《计算机数值方法》 课程设计报告 题目 四阶Runge-Kutta方法 学生姓名 班级 成绩 指导教师 年月日 的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的比较。 三、方法原理及实现 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。 经典的方法是一个四阶的方法，它的计算公式是： 四、计算公式或算法 输入（编写或调用计算的函数文件）， 3．For 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; end y 六、测试数据及结果 用调试好的程序解决如下问题： 应用经典的四阶Runge-Kutta方法解初值问题 取 步骤一：编写函数具体程序. 1.求解解析解程序： dsolve(Dy=(y^2+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(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）步骤二：执行上述Runge-Kutta算法，计算结果为 1.0000 1.5000 2.0000 2.5000 3.0000 -2.0000 -1.4954 -1.3306 -1.2480 -1.1985 -2.0000 -1.5000 -1.3333 -1.2500 -1.2000 0 0.0046 0.0028 0.0020 0.0015 （3）使用Matlab绘图函数“plot(x,y)”绘制问题数值解和解析解的图形。 数值解的图形： plot(x,y) 解析解的图形 plot(x,yy) （4）使用Matlab中的ode45求解，并绘图。 编写函数如下： %ode.m function dy=ode(x,y) dy=(y^2+y)/x; T,Y]=ode45(ode,[1 3],-2); plot(T,Y) 运行结果如下： 结果分析 由图可知此方法与精确解的契合度非常好，基本上与精度解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。 方法改进 同时，由于误差的存在，我们总想尽可能的是误差趋近于零，常用的就是传统的增加取值的个数。最后，我们通过改变步长来进行改进。 具体实现： （1）h=0.1 a=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