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. For each row of XDATA, the specified gridpoint values must be either strictly increasing or strictly decreasing. 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	Description
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 + 2y + 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