资源描述:
《经典Runge-Kutta方法和显式四阶Adams方法》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、微分方程数值解实验报告姓名丁建伟学号200708020211日期2010.10.23实验项目经典Runge-Kutta方法和显式四阶Adams方法指导教师徐强一、上机实验的问题和要求(需求分析):考虑一阶常微分方程初值问题dy/dx=λy,0<=x<=1,y(0)=1,其精确解为y=eλx,求初值问题数值解。要求:(1)λ=-2时,分别以h=1/16,1/32,1/64,1/128为步长时,使用两种方法求解初值问题,给出右端点的误差,分析误差阶。(2)λ=-100时,取h=1/16,用经典Runge-Kutta法求解
2、初值问题,并验证时数值不稳定性。二、程序设计的基本思想,原理和算法描述:算法的基本思想:若已知dy/dx=f(x,y),0<=x<=1;y(0)=1;经典Runge-Kutta法公式:显性四阶Adams公式:算法描述:先设置两个数组y[]和Y[],用来存储结果。经典Runge-Kutta算法:(1)置h=1/n;x=0,y=1;f=my(2)设置一个Runge-Kutta函数,每一步求出对应的x和y值,共循环进行n次,并每次将结果存储到一个数组里y[]。(3)将结果输出,并求在右端点处截断误差。显性四阶Adams:(
3、1)从y0到y3先由经典Runge-Kutta算法求出,并分别存储到Y[0]到Y[3]中;(2)从y4开始,利用前面求出的值,求解y值,并存储到数组相应位置,使用公式为Y[i]=Y[i-1]+(h/24)*(55*m*Y[i-1]-59*m*Y[i-2]+37*m*Y[i-3]-9*m*Y[i-4])。三、主要程序代码或命令:#include#include//因为用到fabs()和exp()函数,加上函数预处理名。#defineM200//申请两个大容量的数组用来存放两种方法每一步
4、求的y值floatm;//将m和h都设为全局变量,后面的分函数也要用到。floath;floatRungeKutta(floaty)//经典Runge-kutta方法的函数,精髓。{floatk1,k2,k3,k4;k1=m*y;k2=m*(y+h*k1/2);k3=m*(y+h*k2/2);k4=m*(y+h*k3);y=y+(h/6)*(k1+2*k2+2*k3+k4);returny;}voidmain()//主函数{intn,i;floaty[M],Y[M];//申请两个数组floaty0,e0,e1;pri
5、ntf("请输入区间等分数n值和参数值m:");scanf("%d,%f",&n,&m);h=1/float(n);y[0]=1.0,Y[0]=1.0;//赋初始值。for(i=1;i<=n;i++)//按照区间等分数,循环执行以下语句n次{if(i<=3)//使用显性4阶Adams方法只能从y4开始求起,故y0到y3都用{//龙塔方法求出。y[i]=RungeKutta(y[i-1]);Y[i]=RungeKutta(Y[i-1]);}else{y[i]=RungeKutta(y[i-1]);Y[i]=Y[i-1]
6、+(h/24)*(55*m*Y[i-1]-59*m*Y[i-2]+37*m*Y[i-3]-9*m*Y[i-4]);}//上公式是显性四阶Adams方法的核心公式。printf("请输出Runge-Kutta算法结果y0:%f和使用显式四阶Adams算法结果y1:%f",y[i],Y[i]);}y0=exp(m*1.0);//求出在区间右端点的精确解e0=fabs(y0-y[n]);//求Runge-Kutta算法在右端点处的截断误差e1=fabs(y0-Y[n]);//求显式四阶Adams算法在右端点处截断误差p
7、rintf("请输出在右端点处的Runge-Kutta算法误差e0=%f和显式四阶Adams算法误差e1=%f",e0,e1);}注意:若想仅使用经典Runge-Kutta算法可修改为for(i=1;i<=n;i++)y[i]=RungeKutta(y[i-1]);//仅用这些就可以了四、调试和运行程序过程中产生的问题及采取的措施:1、编译时出错,若想在主函数和被调用函数都使用一些变量,必须把这些变量设为全局变量。2、编译时,没有注意数据类型转换,如floath;h=1/n;是错误的,因为n是整形的,当n值大于1时,
8、h老为零。应进行模式转换,h=1/float(n);这样才是正确的。3、对浮点数求绝对值时,应使用fabs()函数,而不是abs()。五、运行输出结果及分析:上述程序在VisualC++6.0环境下加以实现。经过多次测试,程序运行正确。例如:分别输入n值和m值:16,-2,运行结果如图所示,图中显示了每一步的值及端点误差。由上图可知:1.在右端