d2r9.cpp

Visual C++ 常用数值算法集 源代码

# include<iostream.h>
# include<math.h>
# include<iomanip.h>

void polint(double xa[],double ya[], int n,double x,double& y,double& dy)
{
	  double den,ho,hp,w,c[10], d[10],dif,dift;
      int i,m,ns = 1;
      dif = fabs(x - xa[1]);
      for (i = 1;i<=n; i++)
	  {
          dift = fabs(x - xa[i]);
          if( dift < dif) 
		  {
              ns = i;
              dif = dift;
		  }
          
          c[i] = ya[i];
          d[i] = ya[i];
      }
      y = ya[ns];
      ns = ns - 1;
      for (m = 1 ;m<=n-1;m++)
	  {
          for (i = 1 ;i<=n-m;i++)
		  {
              ho = xa[i] - x;
              hp = xa[i + m] - x;
              w = c[i + 1] - d[i];
              den = ho - hp;
              if( den == 0.0)
			  {
                 cout<< "pause"<<endl;
                 return;
			  }
              
              den = w / den;
              d[i] = hp * den;
              c[i] = ho * den;
          }
          if (2 * ns < n - m )
              dy = c[ns + 1];
          else
		  {
              dy = d[ns];
              ns = ns - 1;
          }
          y = y + dy;
	  }
}

void polin2(double x1a[],double x2a[],double ya[][6],int m,int n, double x1, double x2,double&y, double& dy)
{
	double yntmp[20], ymtmp[20];
	int j,k;
    for(j = 1; j<=m; j++)
	{
        for(k = 1; k<=n; k++)
            yntmp[k] = ya[j][k];
        polint(x2a, yntmp, n, x2, ymtmp[j], dy);
    }
    polint(x1a, ymtmp, m, x1, y, dy);
}

void main() 
{
    //program d2r9
    //driver for routine polin2
    int n,i,j;
	double pi,x1,x2,f,dy,y,x1a[6],x2a[6],ya[6][6];
    n = 5;
    pi = 3.141593;
    for (i = 1; i<= n; i++)
	{
        x1a[i] = i * pi / n;
        for (j = 1 ;j<= n;j++)
		{
            x2a[j] = 1.0* j / n;
            ya[i][j] = sin(x1a[i]) * exp(x2a[j]);
        }
    }
    //test 2-dimensional interpolation
    cout<<endl;
    cout<<"       x1        x2        f[x]   interpolated  error"<<endl;
    for (i = 1; i<= 4; i++)
	{
        x1 = (-0.1 + i / 5.0) * pi;
        for (j = 1; j<= 4; j++)
		{
            x2 = -0.1 + j / 5.0;
            f = sin(x1) * exp(x2);
            polin2(x1a, x2a, ya, n, n, x1, x2, y, dy);
			cout<<setprecision(6)<<setiosflags(ios::fixed);
            cout<<setw(11)<<x1;
            cout<< setw(11)<<x2;
            cout<< setw(11)<<f;
            cout<< setw(11)<<y;
            cout<< setw(11)<<dy<<endl;
        }
        cout<<"**********************************"<<endl;
    }
}