pgzn.rar_JUT_Simpson求积公式_simpson

#include <iostream.h> #include <math.h> double Y[21],Y1[21],X[21],m[21]; static double u; double first(double x[],double y[],int N); main() { double second(double x[],double y[],int N); double third(double x[],double y[],int N); double forth(int N); int i,j; double x[21],y[21]; cout<<"请您输入插值点数：请输入5或者10或者20:"; int n; cin>>n; switch(n) { case 5: cout<<endl<<" ****** 当n="<<n<<"时: ******"<<endl<<endl; for(j=0;j<=n;j++) X[j]=x[j]=-1+2.0*j/n; first(x,y,n); //调用函数first,计算x[i]所对应的真值，并将真值存于数组Y中； for(j=0;j<=n;j++) Y[j]=y[j]; second(x,y,n); //调用函数second，计算Newton插值多项式的各阶差商，并将其值存于数组Y1中； for(j=0;j<=n;j++) Y1[j]=y[j]; third(x,Y,n); //调用函数third，计算三次样条函数的参数，并将其值存于数组m中； forth(n); //调用函数forth(求x=-1+2*k/100 (k=1~,10,90~,99)时相应的函数值y=f(x),及Nn(x),Sn(x))，E(Nn)，E(Sn); break; case 10: cout<<endl<<" ****** 当n="<<n<<"时: ******"<<endl<<endl; for(j=0;j<=n;j++) X[j]=x[j]=-1+2.0*j/n; first(x,y,n); for(j=0;j<=n;j++) Y[j]=y[j]; second(x,y,n); for(j=0;j<=n;j++) Y1[j]=y[j]; third(x,Y,n); forth(n); break; case 20: cout<<endl<<" ****** 当n="<<n<<"时: ******"<<endl<<endl; for(j=0;j<=n;j++) X[j]=x[j]=-1+2.0*j/n; first(x,y,n); for(j=0;j<=n;j++) Y[j]=y[j]; second(x,y,n); for(j=0;j<=n;j++) Y1[j]=y[j]; third(x,Y,n); forth(n); break; default: cout<<"对不起，请输入5或者10或者20"; break; } return 0; } double first(double x[],double y[],int N) { int j; for(j=0;j<=N;j++) { y[j]=1/(1+25*pow(x[j],2)); cout<<"y["<<j<<"]="<<y[j]<<endl; } return 0; } double second(double x[],double y[],int N) { int i,j; cout<<endl<<"Newton插值多项式的系数："<<endl; for(i=1;i<=N;i++) for(j=N;j>=i;j--) y[j]=(y[j]-y[j-1])/(x[j]-x[j-i]); for(i=0;i<=N;i++) cout<<"y["<<i<<"]="<<y[i]<<endl; return 0; } double third(double x[],double y[],int N) //计算三次样条函数的参数时，还要调用追赶法解三角矩阵的TSS算法； { double tss(double a[],double b[],double c[],double M[],int N); int i,k; double a[21],b[21],M[21],h[21],c[21],z[21]; for(i=0;i<=N;i++) m[i]=y[i]; for(k=1;k<3;k++) for(i=N;i>=k;i--) m[i]=(m[i]-m[i-1])/(x[i]-x[i-k]); h[1]=x[1]-x[0]; for(i=1;i<=N-1;i++) { h[i+1]=x[i+1]-x[i]; c[i]=h[i+1]/(h[i]+h[i+1]); a[i]=1-c[i]; b[i]=2; m[i]=6*m[i+1]; } for(i=0;i<=3;i++) z[i]=y[i]; for(i=1;i<=3;i++) for(k=3;k>=i;k--) z[k]=(z[k]-z[k-1])/(x[k]-x[k-i]); m[0]=-12*h[1]*z[k+1]; for(i=N-3;i<=N;i++) z[i]=y[i]; for(i=N-2;i<=N;i++) for(k=N;k>=i;k--) z[k]=(z[k]-z[k-1])/(x[k]-x[k-i]); m[N]=12*h[N]*z[k+1]; c[0]=-2; b[0]=2; a[N]=-2; b[N]=2; for(i=0;i<=N;i++) M[i]=m[i]; tss(a,b,c,M,N); return 0; } double tss(double a[],double b[],double c[],double M[],int N) { double q[100],y[100],p[100]; int k; q[0]=b[0];y[0]=M[0]; for(k=1;k<=N;k++) { p[k]=a[k]/q[k-1]; q[k]=b[k]-p[k]*c[k-1]; y[k]=M[k]-p[k]*y[k-1]; } m[N]=y[N]/q[N]; cout<<endl<<"三次样条插值的系数："<<endl; cout<<"m["<<N<<"]="<<m[N]<<endl; for(k=N-1;k>0;k--) { m[k]=(y[k]-c[k]*m[k+1])/q[k]; cout<<"m["<<k<<"]="<<m[k]<<endl; } return 0; } double forth(int N) { int i; double Newton(double x,double y,int N,int j); double SN(double x,double y,int N,int j); double x[20],y[20]; for(i=0;i<=19;i++) { if(i<=9) x[i]=-1+2.0*(i+1)/100; if(i>=10) x[i]=-1+2.0*(i+80)/100; } cout<<endl<<endl<<"k=1--10,90--99时："<<endl<<endl; first(x,y,19); //调用函数first，计算出第三小问所给点的函数值； cout<<endl<<"所求牛顿插值的值："<<endl; for(i=0;i<=19;i++) Newton(x[i],y[i],N,i); //调用函数Newton，用Newton法求出(x[i],y[i])点对应的Nn(x)，并求出E(Nn); cout<<endl<<"用Newton法所得最大误差EN="<<u<<endl; cout<<endl<<"所求三次样条插值的值："<<endl; for(i=0;i<=19;i++) SN(x[i],y[i],N,i); //调用函数SN，用三次样条函数法求出(x[i],y[i])电对应的Sn(x))，并计算出E(Sn); cout<<endl<<"用三次样条法所得最大误差ES="<<u<<endl; return 0; } double Newton(double x,double y,int N,int j) { int i; double y1; y1=Y1[N]; for(i=N-1;i>=0;i--) //计算出Nn(x)； y1=Y1[i]+(x-X[i])*y1; cout<<"x="<<x<<"时,y="<<y1<<endl; if(j==0) //先保存第一个E(Nn)在u中,然后将第二个的E(Nn)与u比较，保存大值于u中，再将第三个E(Nn)与u比较，保存大值于u中，以此循环； u=fabs(y-y1); if(fabs(y-y1)>u) u=fabs(y-y1); return 0; } double SN(double x,double y,int N,int j) { int i,k=1; double h,x1,x2,y1; for(i=1;i<=N-1;i++) //找出第三小问所给点的区间； { if(x<X[i]) { k=i; goto next; } else k=i+1; } next: //计算出Sn(x))； h=X[k]-X[k-1]; x1=X[k]-x; x2=x-X[k-1]; y1=(m[k-1]*pow(x1,3)/6+m[k]*pow(x2,3)/6+(Y[k-1]-m[k-1]*pow(h,2)/6)*x1+(Y[k]-m[k]*pow(h,2)/6)*x2)/h; cout<<"x="<<x<<"时，y="<<y1<<endl; if(j==0) //先保存第一个E(Sn)在u中,然后将第二个的E(Sn)与u比较，保存大值于u中，再将第三个E(Sn)与u比较，保存大值于u中，以此循环； u=fabs(y-y1); if(fabs(y-y1)>u) u=fabs(y-y1); return 0; }