《数值分析-常微分数值解的求法C语言》由会员分享，可在线阅读，更多相关《数值分析-常微分数值解的求法C语言（13页珍藏版）》请在金锄头文库上搜索。
1、 本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生学生学号指导教师实验地点实验成绩二 一 六 年 六 月 二 一 六 年 六摘要常微分方程数值解法是关键词：数值解法；常微分 / 1. 实验容和要求常微分方程初值问题有准确解。要求：分别取步长h = 0.1，0.01，0.001，采用改良的Euler方法、4阶经典龙格-库塔RK方法和4阶Adams预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和准确值，并比拟他们的计算准确度。其中多步法需要的初值由经典R-K法提供。运算时，取足以表示计算精度的有效位。2. 算法说明2.1函数与变量说明表1 函数与
2、变量定义函数/变量名类 型说 明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、 外插公式为显式，插公式为隐式。故用插外插公式时
3、需要进展迭代。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
4、*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,d
5、ouble 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+)
6、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
7、=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=
8、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
9、=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、 改良欧拉
10、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
11、);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=10;i+)printf(%0.3f %lf %lf %0.20fn,(a+i)*h,yi,yxi,fabs(yxi-yi);if(selet=5)break;if(selet=6)exit(1);continue;4. 实验结果1、 改良欧拉方法步长x值00.020.040.060.080.10.120.140.160.180.110.952159340.801567920.557746760.257209793-0.04785128-0.301464441-0.455034006-0.47640194-0.355678930.0110.9537942270.8039350080.5599093790.258
