IBM Books

Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

STPINT and DTPINT--Local Polynomial Interpolation

These subroutines perform a polynomial interpolation at specified abscissas, using data points selected from a table of data.

Table 156. Data Types

x, y, t, s, aux Subroutine
Short-precision real STPINT
Long-precision real DTPINT

Syntax

Fortran CALL STPINT | DTPINT (x, y, n, nint, t, s, m, aux, naux)
C and C++ stpint | dtpint (x, y, n, nint, t, s, m, aux, naux);
PL/I CALL STPINT | DTPINT (x, y, n, nint, t, s, m, aux, naux);

On Entry

x
is the vector x of length n, containing the abscissas of the data points used in the interpolations. The elements of x must be distinct and sorted into ascending order. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 156.

y
is the vector y of length n, containing the ordinates of the data points used in the interpolations. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 156.

n
is the number of elements in vectors x and y--that is, the number of data points. Specified as: a fullword integer; n >= 0.

nint
is the number of data points to be used in the interpolation at any given point. Specified as: a fullword integer; 0 <= nint <= n.

t
is the vector t of length m, containing the abscissas at which interpolation is to be done. For optimal performance, t should be sorted into ascending order. Specified as: a one-dimensional array of (at least) length m, containing numbers of the data type indicated in Table 156.

s
See On Return.

m
is the number of elements in vectors t and s--that is, the number of interpolations to be performed. Specified as: a fullword integer; m >= 0.

aux
has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is the storage work area used by this subroutine. Its size is specified by naux.

Specified as: an area of storage, containing numbers of the data type indicated in Table 156. On output, the contents are overwritten.

naux
is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, STPINT and DTPINT dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >= nint+m.

On Return

s
is the vector s of length m, containing the resulting interpolated values; that is, each si is the value of the interpolating polynomial evaluated at ti. Returned as: a one-dimensional array of (at least) length m, containing numbers of the data type indicated in Table 156.

Notes
  1. Vectors x, y, and t must have no common elements with vector s or work area aux; otherwise, results are unpredictable. See Concepts.
  2. The elements of vector x must be distinct and must be sorted into ascending order; that is, x1 < x2 < ...  < xn. Otherwise, results are unpredictable. For details on how to do this, see ISORT, SSORT, and DSORT--Sort the Elements of a Sequence.
  3. The elements of vector t should be sorted into ascending order; that is, t1 <= t2 <= t3 <= ...  <= tm. Otherwise, performance is affected.
  4. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see Using Auxiliary Storage in ESSL.

Function

Polynomial interpolation is performed at specified abscissas, ti for i = 1, m, in vector t, using nint points selected from the following data:

(xj, yj)    for j = 1, n

where:

x1 < x2 < x3 < ... < xn
xj are elements of vector x.
yj are elements of vector y.

The points (xj, yj), used in the interpolation at a given abscissa ti, are chosen as follows, where k = nint/2:

For ti <= xk+1, the first nint points are used.
For ti > xn -nint+k, the last nint points are used.
Otherwise, points h through h+nint-1 are used, where:
xh+k-1 < ti <= xh+k

The interpolated value at each ti is returned in si for i = 1, m. See references [15] and [54]. If n, nint, or m is 0, no computation is performed. For a definition of the polynomial interpolation function performed through a set of data points, see Function.

For STPINT, the Newton divided differences and interpolating values are accumulated in long precision.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

None

Input-Argument Errors
  1. n < 0
  2. nint < 0 or nint > n
  3. m < 0
  4. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value specified in the syntax for this argument. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example shows interpolation using two data points--that is, linear interpolation--at each ti value.

Call Statement and Input
             X   Y   N   NINT   T   S   M   AUX  NAUX
             |   |   |     |    |   |   |    |    |
CALL STPINT( X , Y , 10 ,  2  , T , S , 5 , AUX , 7  )
 
X        =  (0.0, 0.4, 1.0, 1.5, 2.1, 2.6, 3.0, 3.4, 3.9, 4.3)
Y        =  (1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 4.0, 3.0, 2.0, 1.0)
T        =  (-1.0, 0.1, 1.1, 1.2, 3.9)

Output
S        =  (-1.5000, 1.2500, 3.2000, 3.4000, 2.0000)

Example 2

This example shows interpolation using three data points--that is, quadratic interpolation--at each ti value.

Call Statement and Input
             X   Y   N   NINT   T   S   M   AUX  NAUX
             |   |   |     |    |   |   |    |    |
CALL STPINT( X , Y , 10 ,  3  , T , S , 5 , AUX , 8  )
 
X        =  (0.0, 0.4, 1.0, 1.5, 2.1, 2.6, 3.0, 3.4, 3.9, 4.3)
Y        =  (1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 4.0, 3.0, 2.0, 1.0)
T        =  (-1.0, 0.1, 1.1, 1.2, 3.9)

Output
S        =  (-2.6667, 1.2750, 3.2121, 3.4182, 2.0000)


[ Top of Page | Previous Page | Next Page | Table of Contents | Index ]