/*-----------------------------------------------------------------------------*\ | Shepard Method for Bivariate Interpolation shepard2d-demo.cc | | of Scattered Data | | | | Last change: Dec 17, 1997 | | | | Matpack Library Release 1.7.3 | | Copyright (C) 1991-1997 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" //-----------------------------------------------------------------------------// // 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 //-----------------------------------------------------------------------------// double fq (double xx, double yy) { return sqr((xx + 2*yy)/3.0); } double fx (double xx, double yy) { return 2*(xx + 2*yy)/9.0; } 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 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 (int j = 1, k = 0; j <= N; j++) { double yk = double(N-j)/double(N-1); for (int 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 (int 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 (int j = 1; j <= M; j++) { double py = p(j); for (int 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; } //-----------------------------------------------------------------------------//