微分方程数值解上机报告-常微分方程初值问题数值解法.doc

第一章 常微分方程初值问题数值解法 实验一 实验题目 试用(a)欧拉格式 (b)中点格式 (c)预报—校正格式 (d)经典四级四阶R-K格式 编程计算方程： 程序 #includeiostream.h #includestdlib.h #includemath.h const int N=11; double fund(double x,double y); void Euler(double a,double h,double y[]); void Center(double a,double h,double y[]); void YuJiao(double a,double h,double y[]); void SjSj(double a,double h,double y[]); void Adams(double a,double h,double y[]); void main() { double a,b,h,y[N]; int i; char option; cout请输入定义域区间[a,b]：\n; cinab; cout请输入初值y[0]：\n; ciny[0]; h=(b-a)/(N-1); couta:欧拉法 \nb:中心法 \nc:预报-校正格式 \n; coutd:经典四级四阶R-K格式 \n; coute:Adams预报-修正格式 \nf:退出\n; label: cout请选择算法：; cinoption; switch(option) { case a: Euler(a,h,y);break; case b: Center(a,h,y);break; case c: YuJiao(a,h,y);break; case d: YuJiao(a,h,y);break; case e: Adams(a,h,y);break; case f: exit(1);break; default : {cout选择有错,重新选择！\n;goto label;} } cout计算结果为：\n xn\t\tynendl; for(i=0;iN;i++) cout a+i*h\t\ty[i]endl; } void Euler(double a,double h,double y[]) { int i; for(i=1;iN;i++) {y[i]=y[i-1]+h*fund(a+(i-1)*h,y[i-1]);} } void Center(double a,double h,double y[]) { int i; double w; for(i=1;iN;i++) {w=fund(a+(i-1)*h,y[i-1]); y[i]=y[i-1]+h*fund(a+(i-1)*h/2,w*h/2+y[i-1]);} } void YuJiao(double a,double h,double y[]) { int i; double sun; double y_Begin[N]={y[0]}; for(i=1;iN;i++) {y_Begin[i]=y_Begin[i-1]+h*fund(a+(i-1)*h,y_Begin[i-1]);} for(i=1;iN;i++) { while(fabs(y_Begin[i]-sun)0.0001) { sun=y[i-1]+h/2.0*(fund(a+(i-1)*h,y[i-1])+fund(a+i*h,y_Begin[i])); y_Begin[i]=sun; } y[i]=sun; } } void SjSj(double a,double h,double y[]) { int i; double k1,k2,k3,k4; for(i=1;iN;i++) { k1=fund(a+(i-1)*h,y[i-1]); k2=fund(a+(i-1)*h+h/2,y[i-1]+h*k1/2); k3=fund(a+(i-1)*h+h/2,y[i-1]+h*k2/2); k4=fund(a+(i-1)*h+h,y[i-1]+h*k3); y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4); } } void Adams(double a,double h,double y[]) { int i; double me,x[N]; double k1,k