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
