This function evaluates the derivative of a function defined on a rectangular grid using quadratic interpolation.
QD2DR Value of the (IXDER, IYDER) derivative of the function at (X, Y). (Output)
IXDER Order of the x-derivative. (Input)
IYDER Order of the y-derivative. (Input)
X X-coordinate of the point at which the function is to be evaluated. (Input)
Y Y-coordinate of the point at which the function is to be evaluated. (Input)
XDATA Array of
length NXDATA
containing the location of the data points in the
x-direction. (Input)
XDATA must be
increasing.
YDATA Array of
length NYDATA
containing the location of the data points in the
y-direction. (Input)
YDATA must be
increasing.
FDATA Array of
size NXDATA by
NYDATA
containing function values. (Input)
FDATA(I, J) is the value of the
function at (XDATA(I), YDATA(J)).
NXDATA Number
of data points in the x-direction. (Input)
NXDATA must be at
least three.
Default: NXDATA = size (XDATA,1).
NYDATA Number of data
points in the y-direction. (Input)
NYDATA must be at
least three.
Default: NYDATA = size (YDATA,1).
LDF Leading
dimension of FDATA exactly as
specified in the dimension statement of the calling program. (Input)
LDF must be
at least as large as NXDATA.
Default:
LDF = size
(FDATA,1).
CHECK Logical
variable that is .TRUE. if checking of
XDATA and YDATA is required or
.FALSE. if
checking is not required. (Input)
Default: CHECK = .TRUE.
Generic: QD2DR (IXDER, IYDER, X, Y, XDATA, YDATA, FDATA [, ])
Specific: The specific interface names are S_QD2DR and D_QD2DR.
Single: QD2DR (IXDER, IYDER, X, Y, NXDATA, XDATA, NYDATA, YDATA, FDATA, LDF, CHECK)
Double: The double precision fucntion name is DQD2DR.
The function QD2DR
interpolates a table of values, using quadratic polynomials, returning an
approximation to the tabulated function. Let (xi, yj, fij) for i = 1,
, nx and j = 1,
, ny be the tabular data.
Given a point (x, y) at which an interpolated value is desired, we
first find the nearest interior grid point (xi, yj). A bivariate
quadratic interpolant q is then formed using six points near (x,
y). Five of the six points are (xi, yj), (xią1,
yj), and (xi, yją1).
The sixth point is the nearest point to (x, y) of the grid points
(xią1,
yją1).
The value q(p, r)
(x, y) is returned by QD2DR,
where
p = IXDER
and r = IYDER.
1. Informational errors
Type Code
4 6 The XDATA values must be strictly increasing.
4 7 The YDATA values must be strictly increasing.
2. Because quadratic interpolation is used, if the order of any derivative is greater than two, then the returned value is zero.
In this example, the partial derivatives of sin(x + y) at x = y = π/3 are approximated by using QD2DR on a table of size 21 Ũ 42 equally spaced values on the rectangle [0, 2] Ũ [0, 2].
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDF, NXDATA, NYDATA
PARAMETER (NXDATA=21, NYDATA=42, LDF=NXDATA)
!
INTEGER I, IXDER, IYDER, J, NOUT
REAL F, FDATA(LDF,NYDATA), FLOAT, FU, FUNC, PI, Q,&
SIN, X, XDATA(NXDATA), Y, YDATA(NYDATA)
INTRINSIC FLOAT, SIN
EXTERNAL FUNC
! Define function
F(X,Y) = SIN(X+Y)
! Set up X-grid
DO 10 I=1, NXDATA
XDATA(I) = 2.0*(FLOAT(I-1)/FLOAT(NXDATA-1))
10 CONTINUE
! Set up Y-grid
DO 20 I=1, NYDATA
YDATA(I) = 2.0*(FLOAT(I-1)/FLOAT(NYDATA-1))
20 CONTINUE
! Evaluate function on grid
DO 30 I=1, NXDATA
DO 30 J=1, NYDATA
FDATA(I,J) = F(XDATA(I),YDATA(J))
30 CONTINUE
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99998)
! Check XDATA and YDATA
! Get value for PI and set X and Y
PI = CONST('PI')
X = PI/3.0
Y = PI/3.0
! Evaluate and print the function
! and its derivatives at X=PI/3 and
! Y=PI/3.
DO 40 IXDER=0, 1
DO 40 IYDER=0, 1
Q = QD2DR(IXDER,IYDER,X,Y,XDATA,YDATA,FDATA)
FU = FUNC(IXDER,IYDER,X,Y)
WRITE (NOUT,99999) X, Y, IXDER, IYDER, FU, Q, (FU-Q)
40 CONTINUE
!
99998 FORMAT (32X, '(IDX,IDY)', /, 8X, 'X', 8X, 'Y', 3X, 'IDX', 2X,&
'IDY', 3X, 'F (X,Y)', 3X, 'QD2DR', 6X, 'ERROR')
99999 FORMAT (2F9.4, 2I5, 3X, F9.4, 2X, 2F11.4)
END
REAL FUNCTION FUNC (IX, IY, X, Y)
INTEGER IX, IY
REAL X, Y
!
REAL COS, SIN
INTRINSIC COS, SIN
!
IF (IX.EQ.0 .AND. IY.EQ.0) THEN
! Define (0,0) derivative
FUNC = SIN(X+Y)
ELSE IF (IX.EQ.0 .AND. IY.EQ.1) THEN
! Define (0,1) derivative
FUNC = COS(X+Y)
ELSE IF (IX.EQ.1 .AND. IY.EQ.0) THEN
! Define (1,0) derivative
FUNC = COS(X+Y)
ELSE IF (IX.EQ.1 .AND. IY.EQ.1) THEN
! Define (1,1) derivative
FUNC = -SIN(X+Y)
ELSE
FUNC = 0.0
END IF
RETURN
END
(IDX,IDY)
X Y IDX
IDY F (X,Y)
QD2DR ERROR
1.0472
1.0472 0 0
0.8660 0.8661
-0.0001
1.0472 1.0472 0
1 -0.5000
-0.4993 -0.0007
1.0472 1.0472
1 0
-0.5000 -0.4995
-0.0005
1.0472 1.0472 1
1 -0.8660
-0.8634 -0.0026
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |