shepard1.txt

shepard shepard shepard shepard

/*-----------------------------------------------------------------------------*\
| Shepard Method for Bivariate Interpolation		      shepard2d-demo.cpp | 
| of Scattered Data 								|
|                                                                               |
| Last change: Sep 12, 2004							|
|                                                                               |
| Matpack Library Release 1.9.0                                                 |
| Copyright (C) 1991-2004 by Berndt M. Gammel                                   |
|                                                                               |
| Permission to  use, copy, and  distribute  Matpack  in  its  entirety and its | 
| documentation for  non-commercial purpose and  without fee is hereby granted, | 
| provided that this license information and copyright notice appear unmodified | 
| in all copies. This software is provided 'as is'  without  express or implied | 
| warranty.  In no event will the author be held liable for any damages arising | 
| from the use of this software.						|
| Note that distributing Matpack 'bundled' in with any product is considered to | 
| be a 'commercial purpose'.							|
| The software may be modified for your own purposes, but modified versions may | 
| not be distributed without prior consent of the author.			|
|                                                                               |
| Read the  COPYRIGHT and  README files in this distribution about registration	|
| and installation of Matpack.							|
|                                                                               |
\*-----------------------------------------------------------------------------*/

#include "matpack.h"
#include "mpshepard2d.h"
using namespace MATPACK;

//-----------------------------------------------------------------------------//
// This program tests the Shepard scattered data interpolation
// by printing the maximum errors associated
// with interpolated values and gradients on a 10 by 10
// uniform grid in the unit square.  The data set consists
// of 36 nodes with data values taken from a quadratic func-
// tion for which the method is exact.  The ratio of maximum
// interpolation error relative to the machine precision is
// also printed.  This should be o(1).  The interpolated
// values from Shepard_2D_Value and Shepard_2D_Gradient 
// are compared for agreement.
//
// References:
// -----------
//   Robert J. Renka, 
//   Transactions on Mathematical Software, Vol. 14, No. 2, 149 (1988),
//   Algorithm 660, QSHEP2D: Quadratic Shepard Method for Bivariate 
//   Interpolation of Scattered Data.
//
// See also:
// ---------
//   Robert J. Renka,
//   Transactions on Mathematical Software, Vol. 14, No. 2, 139 (1988),
//   Multivariate Interpolation of Large Sets of Scattered Data.
//
//   Robert J. Renka, 
//   Transactions on Mathematical Software, Vol. 14, No. 2, 151 (1988),
//   Algorithm 661, QSHEP3D: Quadratic Shepard Method for Trivariate 
//   Interpolation of Scattered Data.
//-----------------------------------------------------------------------------//

//-----------------------------------------------------------------------------//
// quadratic test function and its partial derivatives
//-----------------------------------------------------------------------------//

static double fq (double xx, double yy) { return sqr((xx + 2*yy)/3.0); }
static double fx (double xx, double yy) { return 2*(xx + 2*yy)/9.0;   }
static double fy (double xx, double yy) { return 4*(xx + 2*yy)/9.0;   }

//-----------------------------------------------------------------------------//

int main (void)
{
  const double eps = 2*DBL_EPSILON;		// precision

  const int N = 6,				// size of test problem
  	    n = N*N, 				// vector dimension
            nq = 13, nw = 19, nr = 3;		// interpolation parameters
  int i,j,k;

  MpShepard2d Interpolator;			// instance of interpolator
  Vector x(1,n), y(1,n), f(1,n); 		// vectors of data points

  // generate a N by N grid of nodes in the unit square
  for (j = 1, k = 0; j <= N; j++) {
    double yk = double(N-j)/double(N-1);
    for (i = 1; i <= N; i++) {
      k++;
      x(k) = double(i-1)/double(N-1);
      y(k) = yk; 
      f(k) = fq(x(k),y(k));
    }
  }
  
  // compute parameters defining the interpolant
  int status;
  status = Interpolator.Interpolate(x,y,f,nq,nw,nr);
  if (status) Matpack.Error("Interpolate: error = %d", status);
 
  // Generate a M by M uniform grid of interpolation points (p(i),p(j)) 
  // in the unit square. The four corners coincide with nodes
  const int M = N*10/6;
  Vector p(1,M);
  for (i = 1; i <= M; i++) p(i) = double(i-1)/double(M-1);
  
  // compute interpolation errors and test for agreement in the
  // q values returned by MpShepard2d::GetValue and MpShepard2d::GetGradient.
  double eq = 0.0, eqx = 0.0, eqy = 0.0, q1, q, qx, qy;
 
  for (j = 1; j <= M; j++) {
    double py = p(j);
    
    for (i = 1; i <= M; i++) {
      double px = p(i);

      // calculate interpolated value at (px,py)
      status = Interpolator.GetValue(px,py,q1);
      if (status) Matpack.Error("GetValue: error = %d", status);
      
      // calculate gradient and interpolated value at (px,py)
      status = Interpolator.GetGradient(px,py,q,qx,qy);
      if (status) Matpack.Error("GetGradient: status = %d", status);
      
      // check size of error - error if values differ by a relative
      // amount greater than 3*eps
      if (abs(q1-q) > 3.0*abs(q)*eps) 
	Matpack.Error("interpolated values q1 (GetValue) "
		      "and q2 (GetGradient) differ\n"
		      "q1 = %.15g\nq2 = %.15g", q1,q);
      
      eq  = max(eq,abs(fq(px,py)-q));
      eqx = max(eqx,abs(fx(px,py)-qx));
      eqy = max(eqy,abs(fy(px,py)-qy));
    }
  }
  
  // print errors and the ratio eq/eps.
  cout << endl
       << "Maximum absolute errors in the interpolant q and partial" << endl
       << "derivatives (qx,qy) relative to machine precision eps" << endl 
       << endl
       << " function     max error    max error/eps" << endl
       << "    q       " << setw(12) << eq  << setw(12) << eq/eps  << endl
       << "    qx      " << setw(12) << eqx << setw(12) << eqx/eps << endl
       << "    qy      " << setw(12) << eqy << setw(12) << eqy/eps << endl
       << endl;
}

//-----------------------------------------------------------------------------//