数值分析-常微分数值解的求法C语言.doc

本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二 一 六 年 六 月 二 一 六 年 六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.关键词：数值解法；常微分121. 实验内容和要求常微分方程初值问题有精确解。要求：分别取步长h = 0.1，0.01，0.001，采用改进的Euler方法、4阶经典龙格-库塔RK方法和4阶Adams预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K法提供。运算时，取足以表示计算精度的有效位。2. 算法说明2.1函数及变量说明表1 函数及变量定义函数/变量名类 型说 明fdouble要求的微分方程p_fdouble微分方程原函数Euler，Euler_ProK_R,Adams double欧拉，改进欧拉，经典K_R，4阶Adams预测-校正方法MainVoid主函数1、 欧拉方法： （i=0，1，2，3，.n-1） (其中a为初值)2、 改进欧拉方法：（i=0，1，2.n-1） (其中a为初值)3、 经典K-R方法： 4、4阶adams预测-校正方法预测：校正： Adsms内插外插公式联合使用称为Adams预测-校正系统，利用外插公式计算预测，用内插公式进行校正。计算时需要注意以下两点：1、 外插公式为显式，内插公式为隐式。故用内插外插公式时需要进行迭代。2、 这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta格式为它提供初始值，依据局部截断误差公式得：进一步将Adams预测-校正系统加工成下列方案：运用上述方案计算时，要先一步的信息和更前一步的信息因此在计算机之前必须给出初值和，可用其他单步法来计算，则一般令他为零。3. 源程序#include#include#include#include/微分方程double f(double x,double y)return (-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x);/return x*pow(y,-2)*2.0/3.0;/原函数double p_f(double x)return x*x*exp(-x)+cos(2*x);/return pow(1.0+pow(x,2),1.0/3.0);/求精确解double* Exact(double a,double b,double h)int i;double c = (b-a)/h;double *y = new double c+1;for(i=0;i=c;i+)yi= p_f(a+i*h);return y;/欧拉方法double* Euler(double a,double b,double h,double y0)int i;double c = (b-a)/h;double *y = new double c+1;y0=y0;for(i=1;i=c;i+)yi= yi-1+ h* f(a+(i-1)*h,yi-1);/printf(%fn,f(a+(i-1)*h,yi-1);return y;/改进欧拉方法double* Euler_Pro(double a,double b,double h,double y0)int i;double c = (b-a)/h ,yb;double *y = new double c+1;y0=y0;for(i=1;i=c;i+)yb=yi-1 +h* f(a+(i-1)*h,yi-1);yi = yi-1 + 0.5*h*( f(a+(i-1)*h,yi-1) + f(a+i*h,yb) );return y;/经典K-R方法double* K_R(double a,double b,double h,double y0)double K1,K2,K3,K4,x;int i;double c = (b-a)/h ;double *y = new double c+1;y0=y0;for(i=1;i=c;i+)x=a+(i-1)*h;K1=f(x,yi-1);K2=f(x+0.5*h,yi-1+0.5*h*K1);K3=f(x+0.5*h,yi-1+0.5*h*K2);K4=f(x+h,yi-1+h*K3);yi=yi-1 + h*( K1 + 2*K2 + 2*K3 + K4)/6;return y;/4阶Adams预测-校正方法double *Adams(double a,double b,double h,double y0)int i;double *y;double count = (b-a)/h;double dy100=0,x4=0;double p_0=0,p_1=0,c=0; double *r = new double *count+1;/动态初始化，储存预测值和校值for(i=0;i=count;i+)ri = new double 2;if(r=NULL)printf(内存分配失败n);exit(0);for(i=0;icount;i+)ri0=0;ri1=0;y=K_R(a,b,h,y0);for(i=0;i4;i+)xi=a+h*i;dyi=f(xi,yi);ri0=yi;i=3;while(x3b)x3 = x3+h;p_1=y3+h*(55*dy3-59*dy2+37*dy1-9*dy0)/24; /预估c=p_1+251*(c-p_0)/270;/修正ri0=c;c=f(x3,c);/求fc=y3+h*(9*c+19*dy3-5*dy2+dy1)/24;/校正y3=c-19*(c-p_1)/270;/修正ri1 = y3;dy0=dy1;dy1=dy2;dy2=dy3;dy3=f(x3,y3);p_0=p_1;i+;return r;void main()int i , selet;double a=0, b=2, h=0.001, y0=1;while(1)printf(请输入下限， 上限， 步长 ， 初值：(用空格隔开)n);scanf(%lf %lf %lf %lf,&a,&b,&h,&y0);double c = (b-a)/h;double *yx,*y;yx=Exact(a,b,h);while(1)printf(1、欧拉 2、 改进欧拉 3、经典K-R4、4阶Adams预测-校正 5、修改数据 6、退出n);scanf(%d,&selet);system(cls);if(selet=1)y=Euler(a,b,h,y0);printf(欧拉方法：n);if(selet=2)y=Euler_Pro(a,b,h,y0);printf(改进欧拉方法：n);if(selet=3)y=K_R(a,b,h,y0);printf(经典K-R方法(4阶龙格-库塔方法)：n);if(selet=4)double *r;r=Adams(a,b,h,y0);printf(4阶Adams预测-校正方法:n);printf(x值 预估值 校正值 误差n);printf(%0.3f %lf - -n,0,r00);printf(%0.3f %lf - -n,0.1,r10);printf(%0.3f %lf - -n,0.2,r20);for(i=3;i=10;i+)printf(%0.3f %f %f %0.20fn,a+i*h,ri0,ri1,fabs(ri0-ri1);if(selet=1|selet=2|selet=3)printf(x值 计算值 准确值 误差n);for(i=0;i改进欧拉欧拉方法4阶adams预测校正方法2， 同一个方法的情况下，步长越小 结果准确度越高；参考文献：1 朱建新，李有法，数值计算方法(第三版)，高等教育出版社；2Adam预测-校正系统/view/56260294daef5ef7ba0d3cd9学