interpolate.cpp

能够进行三点插值、牛顿插值和Hermit插值的小程序。

///////////////////////////////////////////////////
// This program is for interpolation.
////////////////////////////////////////////////////
# include <iostream.h>
# include <stdio.h>
# include <math.h>

/////////////////////////////////////////////////////////////////////////////
double ThreePointsTnterpolate(double px[],double py[],int n, double x);
double NewtonInterpolate(double px[], double py[], int n, double x);
double HermiteInterInterpolate(double px[], double py[], double dpy[],int n, double x);
/////////////////////////////////////////////////////////////////////////////

void main()
 { 
    double  t=0.356,z1=0.0,z2=0.0,z3=0.0;
	double  x[10]={0.1,0.15,0.30,0.45,0.55,0.60,0.70,0.85,0.90,1.00};
	double  y[10]={0.904937,0.860708,0.740818,0.637628,0.576950,0.548812,
		           0.496585,0.427415,0.406570,0.367879};
	double dy[10]={-0.904937,-0.860708,-0.740818,-0.637628,-0.576950,
		           -0.548812,-0.496585,-0.427415,-0.406570,-0.367879};


	 for (int i=0; i<5; i++)
    {
        printf("px[%d]=%f ; py[%d]=%f ;&px[i]=%d ;&py[i]=%d \n",i,x[i],i,y[i],&x[i],&y[i]);
    }
	printf("please input zhe point t: \n");
	//cin>>t;
	z1=ThreePointsTnterpolate(x,y,5,t);
	printf("x1=%6.3f,      f1(x)=%f\n",t,z1);
	cout<<"\n";
	z2=NewtonInterpolate(x, y, 5, t);
	printf("x2=%6.3f,      f2(x)=%f\n",t,z2);
	cout<<"\n";
	z3=HermiteInterInterpolate(x,y,dy,5, t);
    printf("x3=%6.3f,      f3(x)=%f\n",t,z3);
	cout<<"\n";
	
}



 double ThreePointsTnterpolate(double px[],double py[],int n, double x)
 {
	 //px,py;插值点(Xi,Yi) n:插值点个数  x:待计算的点
	 
	 int i,j,k,m;
	 double z,s;
	 z=0.0;
	 if(n<1)
	 {
		 return(z);
	 }
	 if(n==1)
	 {
		 z=py[0];return (z);
	 }
	 if(n==2)
	 {
		 if(px[0]!=px[1])
		 {
			 z=(py[0]*(x-px[1])-py[1]*(x-px[0]))/(px[0]-px[1]);
			 return(z);
		 }
		 else
			 printf("The two x value are equal. \n");
	 }
	 else 
	 {
		 if (x<=px[1])
		 {
			 k=0;
			 m=2;
		 }
		 else if(x>=px[n-2])
		 {
			 k=n-3;m=n-1;
		 }
		 else 
		 {
		     k=1;m=n;
			 while ( (m-k)!=1 )
				  {
					  i=(k+m)/2;
				      if( x<px[i-1] )
					       m=i;
				      else 
					       k=i;
				  }
			  k=k-1;
			  m=m-1;
			  if ( fabs(x-px[k])< fabs(x-px[m]) ) 
				  k=k-1;
			  else
				  m=m+1;
		 }
		 z=0.0;
		 for(i=k;i<=m;i++)
		 {
			 s=1.0;
			 for(j=k;j<=m;j++)
			 
				 if(j!=i) 
					 s=s*(x-px[j])/(px[i]-px[j]);
			 z=z+s*py[i];
			
		 }
	}
	 return(z);
}

double  NewtonInterpolate(double px[], double py[], int n, double x)
{   //px,py;插值点(Xi,Yi) n:插值点个数  x:待计算的点
    double *g = new double [n];

    //计算差商表个g[k],k=0,1,...,n
	
    for (int i=0; i<=n-1; i++)
    {
        g[i] = py[i];
		 
	}

	for (int k=1; k<=n-1; k++)
    {
	    for ( i=n-1; i>=k; i--)
        {
            g[i] = (g[i] - g[i-1]) / (px[i] - px[i-k]);
			// g[i]用来暂时存放f[X(i-k), X(i)]
		}
    }
    
    double z = 1;//基函数
    double newton = py[0];
    for (i=1; i<n; i++)
    {
        z = (x - px[i-1]) * z;//基函数的递推式
        newton = newton +  z * g[i];
    }

    delete [] g;
    return newton;
}

double HermiteInterInterpolate(double px[], double py[], double dpy[],int n, double x)
{
	int i=0,j=0;
	double z=0.0,p=0.0,q=0.0,s=0.0;
	
	for (i=0;i<n;i++)
	{
		s=1.0;
		for (j=0;j<n;j++)
			if(i!=j)
				   s=s*(x-px[j])/(px[i]-px[j]);
		s=s*s;
		p=0.0;
		for (j=0;j<n;j++)
			if(i!=j)
				p=p+1.0/(px[i]-px[j]);
		q=py[i]+(x-px[i])*(dpy[i]-2.0*py[i]*p);
		z=z+q*s;
	}
    return z;
}