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

SPINT and DPINT--Polynomial Interpolation

These subroutines compute the Newton divided difference coefficients and perform a polynomial interpolation through a set of data points at specified abscissas.

Table 155. Data Types

Syntax

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. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 155.

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 155.

n

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

c

is the vector c of length n, where:

If ninit <= 0, all elements of c are undefined on entry.

If ninit > 0, c contains the Newton divided difference coefficients, c_j for j = 1, ninit, for the interpolating polynomial through the data points (x_j,y_j) for j = 1, ninit. If ninit < n, the values of c_j for j = ninit+1, n are undefined.

Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 155.

ninit

indicates the following:

If ninit <= 0, this is the first call to this subroutine with the data in x and y; therefore, none of the Newton divided difference coefficients in c have been initialized.

If ninit > 0, a previous call to this subroutine was made with the data points (x_j, y_j) for j = 1, ninit, where:

• If ninit = n, all the Newton divided difference coefficients in c were computed for the data points. No additional coefficients are computed on this entry.
• If ninit < n, the first ninit Newton divided difference coefficients in c were computed for the data points (x_j, y_j) for j = 1, ninit. The coefficients are updated for the additional data points (x_j, y_j) for j = ninit+1, n on this entry.

Specified as: a fullword integer; ninit <= n.

t

is the vector t of length m, containing the abscissas at which interpolation is to be done. Specified as: a one-dimensional array of (at least) length m, containing numbers of the data type indicated in Table 155.

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.

On Return

c

is the vector c of length n, containing the coefficients of the Newton divided difference form of the interpolating polynomial through the data points (x_j,y_j) for j = 1, n. Returned as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 155.

ninit

is the number of coefficients, n, in output vector c. (If you call this subroutine again with the same data, this value should be specified for ninit.) Returned as: a fullword integer; ninit = n.

s

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

Notes

1. In your C program, argument ninit must be passed by reference.
2. Vectors x, y, and t must have no common elements with vectors c and s, and vector c must have no common element with vector s; otherwise, results are unpredictable.
3. The elements of vector x must be distinct; that is, x_i <> x_j if i <> j for i, j = 1, n.

Function

Polynomial interpolation is performed at specified abscissas, t_i for i = 1, m, in vector t, using the method of Newton divided differences through the data points:

(x_j, y_j) for j = 1, n

where:

x_j are elements of vector x.

y_j are elements of vector y.

The interpolated value at each t_i is returned in s_i for i = 1, m. See references [15] and [54]. The interpolating values returned in s are computed using the Newton divided difference coefficients, as defined in the following section.

The divided difference coefficients, c_j for j = 1, n, are returned in vector c. These coefficients can then be reused on subsequent calls to this subroutine, using the same data points (x_j, y_j), but with new values of t_i. If the number of data points is increased from one call this subroutine to the next, the new coefficients are computed, and the existing coefficients are updated (not recomputed). This feature can be used to test for the convergence of the interpolations through a sequence of an increasingly larger set of points.

The values specified for ninit and m indicate which combination of functions are performed by this subroutine: computing the coefficients, performing the interpolation, or both. If m = 0, only the divided difference coefficients are computed. No interpolation is performed. If n = 0, no computation or interpolation is performed.

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

Newton Divided Differences and Interpolating Values

The Newton divided differences of the following data points:

(x_j, y_j)    for j = 1, n

where x_j <> x_l if j <> l    for j, l = 1, n

are denoted by delta_ky_j for k = 0, 1, 2, ..., n-1 and j = 1, 2, ..., n-k, and are defined as follows:

For k = 0 and 1:

delta₀y_j = y_j    for j = 1, 2, ..., n

delta₁y_j = (y_j+1 - y_j) / (x_j+1 - x_j)    for j = 1, 2, ..., n-1

For k = 2, 3, ..., n-1:

delta_ky_j = (delta_k-1 y_j+1 - delta_k-1y_j) / (x_j+k - x_j)    for j = 1, 2, ..., n-k

The value s of the Newton divided difference form of the interpolating polynomial evaluated at an abscissa t is given by:

s = y_n + (t-x_n) delta₁y_n-1

+ (t-x_n-1) (t-x_n) delta₂y_n-2

+ ...+(t-x₂) (t-x₃) ... (t-x_n) delta_n-1y₁

Therefore, on output, the coefficients in vector c are as follows:

c_n = y_n

c_n-1 = delta₁y_n-1

c_n-2 = delta₂y_n-2

.

.

.

c₁ = delta_n-1y₁

Also, the interpolating values in s, in terms of c, are as follows for i = 1, m:

s_i = c_n + (t_i-x_n) c_n-1

+ (t_i-x_n-1) (t_i-x_n) c_n-2

+ ...

+ (t_i-x₂) (t_i-x₃) ... (t_i-x_n) c₁

Error Conditions

Computational Errors

None

Input-Argument Errors

1. n < 0
2. ninit > n
3. m < 0

Example 1

This example shows a quadratic polynomial interpolation on the initial call with the specified data points; that is, NINIT = 0, and C contains all undefined values. On output, NINIT and C are updated with new values.

Call Statement and Input

            X   Y   N   C  NINIT  T   S   M
            |   |   |   |    |    |   |   |
CALL SPINT( X , Y , 3 , C ,  0  , T , S , 2 )

X        =  (-0.50, 0.00, 1.00)
Y        =  (0.25, 0.00, 1.00)
C        =  ( . , . , . )
T        =  (-0.2, 0.2)

Output

C        =  (1.00, 1.00, 1.00)
NINIT    =  3
S        =  (0.04, 0.04)

Example 2

This example shows a quadratic polynomial interpolation on a subsequent call with the same data points specified in Example 1, but using a different set of abscissas in T. In this case, NINIT = N = 3, and C contains the values defined on output in Example 1. On output here, the values in NINIT and C are unchanged.

Call Statement and Input

            X   Y   N   C  NINIT  T   S   M
            |   |   |   |    |    |   |   |
CALL SPINT( X , Y , 3 , C ,  3  , T , S , 2 )

X        =  (-0.50, 0.00, 1.00)
Y        =  (0.25, 0.00, 1.00)
C        =  (1.00, 1.00, 1.00)
T        =  (-0.10, 0.10)

Output

C        =  (1.00, 1.00, 1.00)
NINIT    =  3
S        =  (0.01, 0.01)

Example 3

This example is the same as Example 2 except that it specifies additional data points on the subsequent call to the subroutine. In this case, 0 < NINIT < N. On output here, the values in NINIT and C are updated. The interpolating polynomial is a degree of 4.

Call Statement and Input

            X   Y   N   C  NINIT  T   S   M
            |   |   |   |    |    |   |   |
CALL SPINT( X , Y , 5 , C ,  3  , T , S , 2 )

X        =  (-0.50, 0.00, 1.00, -1.00, 0.50)
Y        =  (0.25, 0.00, 1.00, 1.10, 0.26)
C        =  (1.00, 1.00, 1.00, . , . )
T        =  (-0.10, 0.10)

Output

C        =  (0.04, -0.06, 1.02, -0.56, 0.26)
NINIT    =  5
S        =  (0.0072, 0.0130)