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

本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六1摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.关键词：数值解法；常微分21.实验内容和要求常微分方程初值问题cos22sin2202(0)1xyyxxxexy有精确解2()cos(2)xyxxex。要求：分别取步长h=0.1，0.01，0.001，采用改进的Euler方法、4阶经典龙格-库塔R－K方法和4阶Adams预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K法提供。运算时，取足以表示计算精度的有效位。2.算法说明2.1函数及变量说明表1函数及变量定义函数/变量名类型说明fdouble要求的微分方程p_fdouble微分方程原函数Euler，Euler_ProK_R,Adamsdouble欧拉，改进欧拉，经典K_R，4阶Adams预测-校正方法MainVoid主函数31、欧拉方法：1()()(,())iiiiyxyxhfxyx（i=0，1，2，3，......n-1）(0)ya(其中a为初值)2、改进欧拉方法：~1~111()()(,())()()[(,())(,())]2(0)iiiiiiiiiiyxyxhfxyxhyxyxfxyxfxyxya（i=0，1，2......n-1）(其中a为初值)3、经典K-R方法：11213243()6(,)(,)22(,)22(,)iiiiiiiiiihyyKfxyhhKfxyKhhKfxyKKfxhyhK4、4阶adams预测-校正方法预测：校正：Adsms内插外插公式联合使用称为Adams预测-校正系统，利用外插公式计算预测，用内插公式进行校正。计算时需要注意以下两点：41、外插公式为显式，内插公式为隐式。故用内插外插公式时需要进行迭代。2、这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息1-nf，2-nf3-nf，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta格式为它提供初始值321,,yyy，依据局部截断误差公式得：进一步将Adams预测-校正系统加工成下列方案：运用上述方案计算1ny时，要先一步的信息nnncpy,,y'n和更前一步的信息1-ny因此在计算机之前必须给出初值1y和11c-p，1y可用其他单步法来计算，11c-p则一般令他为零。3.源程序5#includestdio.h#includestdlib.h#includewindows.h#includemath.h//微分方程doublef(doublex,doubley){return(-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x));//returnx*pow(y,-2)*2.0/3.0;}//原函数doublep_f(doublex){returnx*x*exp(-x)+cos(2*x);//returnpow(1.0+pow(x,2),1.0/3.0);}//求精确解double*Exact(doublea,doubleb,doubleh){inti;doublec=(b-a)/h;double*y=newdouble[c+1];for(i=0;i=c;i++)y[i]=p_f(a+i*h);returny;}//欧拉方法double*Euler(doublea,doubleb,doubleh,doubley0){inti;doublec=(b-a)/h;double*y=newdouble[c+1];y[0]=y0;for(i=1;i=c;i++){y[i]=y[i-1]+h*f(a+(i-1)*h,y[i-1]);//printf(%f\n,f(a+(i-1)*h,y[i-1]));}returny;}//改进欧拉方法double*Euler_Pro(doublea,doubleb,doubleh,doubley0){inti;6doublec=(b-a)/h,yb;double*y=newdouble[c+1];y[0]=y0;for(i=1;i=c;i++){yb=y[i-1]+h*f(a+(i-1)*h,y[i-1]);y[i]=y[i-1]+0.5*h*(f(a+(i-1)*h,y[i-1])+f(a+i*h,yb));}returny;}//经典K-R方法double*K_R(doublea,doubleb,doubleh,doubley0){doubleK1,K2,K3,K4,x;inti;doublec=(b-a)/h;double*y=newdouble[c+1];y[0]=y0;for(i=1;i=c;i++){x=a+(i-1)*h;K1=f(x,y[i-1]);K2=f(x+0.5*h,y[i-1]+0.5*h*K1);K3=f(x+0.5*h,y[i-1]+0.5*h*K2);K4=f(x+h,y[i-1]+h*K3);y[i]=y[i-1]+h*(K1+2*K2+2*K3+K4)/6;}returny;}//4阶Adams预测-校正方法double**Adams(doublea,doubleb,doubleh,doubley0){inti;double*y;doublecount=(b-a)/h;doubledy[100]={0},x[4]={0};doublep_0=0,p_1=0,c=0;double**r=newdouble*[count+1];//动态初始化，储存预测值和校值for(i=0;i=count;i++)r[i]=newdouble[2];if(r==NULL){printf(内存分配失败\n);exit(0);7}for(i=0;icount;i++){r[i][0]=0;r[i][1]=0;}y=K_R(a,b,h,y0);for(i=0;i4;i++){x[i]=a+h*i;dy[i]=f(x[i],y[i]);r[i][0]=y[i];}i=3;while(x[3]b){x[3]=x[3]+h;p_1=y[3]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24;//预估c=p_1+251*(c-p_0)/270;//修正r[i][0]=c;c=f(x[3],c);//求fc=y[3]+h*(9*c+19*dy[3]-5*dy[2]+dy[1])/24;//校正y[3]=c-19*(c-p_1)/270;//修正r[i][1]=y[3];dy[0]=dy[1];dy[1]=dy[2];dy[2]=dy[3];dy[3]=f(x[3],y[3]);p_0=p_1;i++;}returnr;}voidmain(){inti,selet;doublea=0,b=2,h=0.001,y0=1;while(1){printf(请输入下限，上限，步长，初值：(用空格隔开)\n);scanf(%lf%lf%lf%lf,&a,&b,&h,&y0);doublec=(b-a)/h;double*yx,*y;yx=Exact(a,b,h);8while(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,r[0][0]);printf(%0.3f%lf--\n,0.1,r[1][0]);printf(%0.3f%lf--\n,0.2,r[2][0]);for(i=3;i=10;i++)printf(%0.3f%f%f%0.20f\n,a+i*h,r[i][0],r[i][1],fabs(r[i][0]-r[i][1]));}if(selet==1||selet==2||selet==3){printf(x值计算值准确值误差\n);for(i=0;i=10;i++)9printf(%0.3f%lf%lf%0.20f\n,(a+i)*h,y[i],yx[i],fabs(yx[i]-y[i]));}if(selet==5)break;if(selet==6)exit(1);continue;}}}4.实验结果1、改进欧拉方法2、经典K-R方法步长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.258360654-0.048262298-0.303651576-0.458854025-0.481390334-0.3611365130.00110.9538100650.8039576890.5599297390.258370912-0.048267344-0.303673825-0.458891909-0.481439197-0.361189485精确解10.9538102240.8039579170.5599299440.258371015-0.048267395-0.30367405-0.458892291-0.48143969-0.361190019步长x值00.020.040.060.080.10.120.140.160.180.110.9538102070.803957910.5599300270.258371272-0.048266908-0.303673326-0.458891378-0.481438681-0.361189040.0110.9538102240.8039579170.5599299440.258371015-0.048267395-0.30367405-0.458892291-0.48143969-0.3611900180.00110.9538102240.8039579170.5599299440.258371015-0.048267395-0.30367405-0.458892291-0.48143969-0.361190019精确解10.9538102240.8039579170.5599299440.258371015-0.048267395-0.30367405-0.458892291-0.48143969-0.361190019103、4阶adams预测-校正法5.对比分析欧拉方法改进欧拉方法经典K-R方法4阶adams预测-校正方法0.10.0380189284290141510.0025896560125808660.000