Chapter 3: Interpolation and Approximation

SURFND

Performs multidimensional interpolation and differentiation for up to 7 dimensions.

The dimension, n, of the problem is determined by the rank of FDATA, and cannot be greater than seven. The number of gridpoints in the i-th direction, di, is determined by the corresponding dimension for FDATA.

Function Return Value

SURFND — Interpolated value of the function.

Required Arguments

X — Array of length n containing the point at which interpolation is to be done.   (Input)
An interpolant is to be calculated at the point:

(X1, X2, …, Xn)

XDATA — Array of size n by max(d1, …, dn) giving the gridpoint values for the function to be interpolated.   (Input)
The gridpoints need not be uniformly spaced.   See FDATA for more details.

FDATA  —   n dimensional array, dimensioned  d1× d2× …× dn giving the values at the gridpoints of the function to be interpolated.  (Input) 
FDATA(i, j, k, …) is the value of the function at

(XDATA1,i, XDATA2,j, XDATA3,k, …)

for i =1, …, d1,  j =1, …, d2, k=1, …, d3,  …

Optional Arguments

NDEG   — Array of length n, giving the degree of polynomial interpolation to be used in each dimension.   (Input)
NDEG(i) must be less than or equal to 15.
Default:  NDEG(i) = 5, for i = 1, ,  n.

NDERS  — Maximum order of derivatives to be computed with respect to each variable.   (Input)
NDERS cannot be larger than max (7- n, 2). See DERIV for more details.
Default: NDERS = 0.

DERIV  — n dimensional array, dimensioned (NDERS+1) × (NDERS+1) ×… containing derivative estimates at the interpolation point.  (Output)
DERIV (i,  j, …) will hold an estimate of the derivative with respect to X1 i times, and with respect to X2 j times, etc.  where i = 0, …, NDERS,  j = 0, …, NDERS, ….  The 0-th derivative means the function value, thus, DERIV (0, 0, …) = SURFND.

ERROR  — Estimate of the error in SURFND.   (Output)

FORTRAN 90 Interface

Generic:                              SURFND (X,XDATA,FDATA [,…])

Specific:                             The specific interface names are Sn_SURFND and Dn_SURFND, where “n” indicates the dimension of the problem (n = 1, 2, 3, 4, 5, 6 or 7).

Description

The function SURFND interpolates a function of up to 7 variables, defined on a (possibly nonuniform) grid.  It fits a polynomial of up to degree 15 in each variable through the grid points nearest the interpolation point. Approximations of partial derivatives are calculated, if requested. If derivatives are desired, high precision is strongly recommended. For more details, see Krogh (1970).

Comments

Informational errors

Type Code

3         1                  NDERS is too large, it has been reset to  max(7- n,2).

3         2                  The interpolation point is outside the domain of the table, so extrapolation is used.

4         3                  Too many derivatives requested for the polynomial degree used.

4         4                  One of the polynomial degrees requested is too large for the number of gridlines in that direction.

Example

The 3D function f(x, y, z) = exp(x + 2 y+ 3z), defined on a 20 by 30 by 40 uniform grid, is interpolated.

 

      USE SURFND_INT

      USE UMACH_INT

      IMPLICIT NONE

     

      INTEGER, PARAMETER :: N=3, ND1=20, ND2=30, ND3=40, NDERS=1

      REAL         X(N),DEROUT(0:NDERS,0:NDERS,0:NDERS), &

                   XDATA(N,MAX(ND1,ND2,ND3)),FDATA(ND1,ND2,ND3), &

                   ERROR,XX,YY,ZZ,TRUE,RELERR,YOUT

      INTEGER      NDEG(N),I,J,K,NOUT

      CHARACTER*1  ORDER(3)

      INTRINSIC    EXP, MAX

     

!                             20 by 30 by 40 uniform grid used for

!                             interpolation of F(x,y,z) = exp(x+2*y+3*z)

      NDEG(1) =  8

      NDEG(2) =  7

      NDEG(3) =  9

     

      DO I=1,ND1

         XDATA(1,I) = 0.05*(I-1)

      END DO

     

      DO J=1,ND2

         XDATA(2,J) = 0.03*(J-1)

      END DO

     

      DO K=1,ND3

         XDATA(3,K) = 0.025*(K-1)

      END DO

     

      DO I=1,ND1

         DO J=1,ND2

            DO K=1,ND3

               XX = XDATA(1,I)

               YY = XDATA(2,J)

               ZZ = XDATA(3,K)

               FDATA(I,J,K) =  EXP(XX+2*YY+3*ZZ)

            END DO

         END DO

      END DO

     

!                             Interpolate at (0.18,0.43,0.35)

      X(1) =  0.18

      X(2) =  0.43

      X(3) =  0.35

     

!                             Call SURFND

      YOUT = SURFND(X,XDATA,FDATA,NDEG=NDEG,DERIV=DEROUT,ERROR=ERROR, &

        NDERS=NDERS)

       

!                             Output F,Fx,Fy,Fz,Fxy,Fxz,Fyz,Fxyz at

!                             interpolation point

      XX = X(1)

      YY = X(2)

      ZZ = X(3)

      CALL UMACH (2, NOUT)

      WRITE(NOUT, 10) YOUT,ERROR

              

      DO K=0,NDERS

         DO J=0,NDERS

            DO I=0,NDERS

               ORDER(1:3) = ' '

               IF (I.EQ.1) ORDER(1) = 'x'     

               IF (J.EQ.1) ORDER(2) = 'y'

               IF (K.EQ.1) ORDER(3) = 'z'

               TRUE = 2**J*3**K*EXP(XX+2*YY+3*ZZ)

               RELERR = (DEROUT(I,J,K)-TRUE)/TRUE

               WRITE(NOUT, 20) ORDER,DEROUT(I,J,K),TRUE,RELERR

            END DO

         END DO

      END DO

   10 FORMAT (' EST. VALUE = ',F10.6,', EST. ERROR = ',E11.3,//, &

              11X,'Computed Der.',5X,'True Der.',4X,'Rel. Err')

   20 FORMAT (2X,'F',3A1,2F15.6,E15.3)

      END

Output

 

EST. VALUE =   8.084915, EST. ERROR =   0.419E-05

 

           Computed Der.     True Der.    Rel. Err

  F          8.084915       8.084914      0.118E-06

  Fx         8.084907       8.084914     -0.944E-06

  F y       16.169882      16.169828      0.330E-05

  Fxy       16.171101      16.169828      0.787E-04

  F  z      24.254705      24.254742     -0.149E-05

  Fx z      24.255133      24.254742      0.161E-04

  F yz      48.505203      48.509483     -0.882E-04

  Fxyz      48.464718      48.509483     -0.923E-03



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260