c++插值代码

//数值计算实验插值#includeiostream#includecmath#includecstdlib#includegsl/gsl_errno.h#includegsl/gsl_spline.h#includemylib/MyClock.h#includemylib/InitData.husingnamespacestd;//要插值的函数,为sin(x)，可改。doublef(doublex){returnsin(x);}//产生函数区间为[a,b]上的n个点的插值函数数值，x,y分别存入横纵座标boolInitFun(constdouble&a,constdouble&b,constint&n,vectordouble&x,vectordouble&y){//如果a大于b返回错误。if(ab){returnfalse;}//计算x值间隔大小，将[a,b]均分为n-1份，这样连端点，一共n个点doublegap=(b-a)/(n-1.0);x.resize(n);y.resize(n);for(inti=0;in;i++){x[i]=a+gap*i;y[i]=f(x[i]);}returntrue;}//根据x,y,n，计算n阶拉格朗日插值y1=f(x1)boolLagrange(constvectordouble&x,constvectordouble&y,constdouble&x1,double&y1){y1=0.0;intn=x.size();for(intj=0;jn;j++){doubleL=1.0;for(inti=0;in;i++){if(i!=j){L*=(x1-x[i])/(x[j]-x[i]);}}y1+=L*y[j];}returntrue;}//根据x,y,n，计算n阶牛顿插值y1=f(x1),//参考(constvectordouble&x,constvectordouble&y,constdouble&x1,double&y1){intn=x.size();//计算差商vectordoublem_diff(n);//先存放0阶差商for(inti=0;in;i++){m_diff[i]=y[i];}for(inti=0;in-1;i++){for(intj=n-1;ji;j--){m_diff[j]=(m_diff[j]-m_diff[j-1])/(x[j]-x[j-1-i]);}}doubletemp=1.0;y1=m_diff[0];for(inti=0;in-1;i++){temp=temp*(x1-x[i]);y1+=temp*m_diff[i+1];}returntrue;}//根据x,y,n，调用gsl函数计算n阶插值y1=f(x1),//参考gsl文档boolGSLInterp(vectordouble&x,vectordouble&y,double&x1,double&y1){intn=x.size();double*gx,*gy;try{gx=newdouble[n];gy=newdouble[n];}catch(bad_alloc){cout内存分配失败!endl;returnfalse;}for(inti=0;in;i++){gx[i]=x[i];gy[i]=y[i];}gsl_interp_accel*acc=gsl_interp_accel_alloc();gsl_spline*spline=gsl_spline_alloc(gsl_interp_cspline,n);gsl_spline_init(spline,gx,gy,n);y1=gsl_spline_eval(spline,x1,acc);gsl_spline_free(spline);gsl_interp_accel_free(acc);delete[]gx;delete[]gy;returntrue;}intmain(){intn;doublea,b;cout请输入插值区间:;cinab;cout请输入样本量:;cinn;if(ab||n=0){cout输入错误!endl;exit(1);}vectordoublex(n);vectordoubley(n);//产生插值数据if(!InitFun(a,b,n,x,y)){cout产生插值数据出错!endl;exit(1);}/*//调试InitFun函数用的else{cout产生的插值数据为:endl;cout.precision(10);cout.width(12);coutx=endl;ShowData(x,cout);couty=endl;ShowData(y,cout);}*///设定输出精度cout.precision(8);cout.width(10);//测试插值函数/*doublex1,y1;cout输入插值区间[a,b]内要插值的x值:;cinx1;if(x1a||x1b){cout输入错误!endl;exit(1);}//if(Lagrange(x,y,x1,y1))//if(Newton(x,y,x1,y1))if(GSLInterp(x,y,x1,y1)){//cout拉格朗日插值结果为:y1endl;//cout牛顿插值结果为:y1endl;coutgsl函数插值结果为:y1endl;//cout所求的精确值为:sqrt(abs(x1))endl;cout所求的精确值为:sin(x1)endl;}else{cout插值失败!endl;}*///批量插值intm;cout输入要插值的数据数量:;cinm;if(m=0){cout输入数值必须大于0endl;exit(1);}vectordoublex2(m);vectordoubley2(m);//拉格朗日插值结果vectordoubley3(m);//牛顿插值结果vectordoubley4(m);//gsl插值结果doublegap=(b-a)/(m-1);for(inti=0;im;i++){x2[i]=a+gap*i;if(Lagrange(x,y,x2[i],y2[i])!=true){cout拉格朗日插值失败，程序退出!endl;exit(1);}if(Newton(x,y,x2[i],y3[i])!=true){cout牛顿插值失败，程序退出!endl;exit(1);}if(GSLInterp(x,y,x2[i],y4[i])!=true){coutgsl插值失败，程序退出!endl;exit(1);}}cout插值情况为:endl;coutx=endl;ShowData(x2,cout);cout拉格朗日插值结果:endl;couty=endl;ShowData(y2,cout);cout牛顿插值结果:endl;couty=endl;ShowData(y3,cout);coutgsl插值结果:endl;couty=endl;ShowData(y4,cout);return0;}