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

SCSINT and DCSINT--Cubic Spline Interpolation

These subroutines compute the coefficients of the cubic spline through a set of data points and evaluate the spline at specified abscissas.

Table 157. Data Types

Syntax

On Entry

x

is the vector x of length n, containing the abscissas of the data points that define the spline. 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 157.

y

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

c

is the matrix C with elements c_jk for j = 1, n and k = 1, 4 that contain the following:

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

If init = 1, c₁₁ contains the spline derivative at x₁.

If init = 2, c₂₁ contains the spline derivative at x_n.

If init = 3, c₁₁ contains the spline derivative at x₁, and c₂₁ contains the spline derivative at x_n.

If init > 3, c contains the coefficients of the spline computed for the data points (x_j,y_j) for j = 1, n on a previous call to this subroutine.

Specified as: an n by (at least) 4 array, containing numbers of the data type indicated in Table 157.

n

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

init

indicates the following, where in those cases for uninitialized coefficients, this is the first call to this subroutine with the data in x and y:

If init <= 0, the coefficients are uninitialized. The second derivatives of the spline at x₁ and x_n are set to zero. (These are free end conditions, also called natural boundary conditions.)

If init = 1, the coefficients are uninitialized. The value in c₁₁ is used as the spline derivative at x₁.

If init = 2, the coefficients are uninitialized. The value in c₂₁ is used as the spline derivative at x_n.

If init = 3, the coefficients are uninitialized. The value in c₁₁ is used as the spline derivative at x₁ and the value in c₂₁ is used as the spline derivative at x_n.

If init > 3, the coefficients in c were computed for data points (x_j, y_j) for j = 1, n on a previous call to this subroutine.

Specified as: a fullword integer. It can have any value.

t

is the vector t of length m, containing the abscissas at which the spline is evaluated. Specified as: a one-dimensional array of (at least) length m, containing numbers of the data type indicated in Table 157.

s

See On Return.

m

is the number of elements in vectors t and s--that is, the number of points at which the spline interpolation is evaluated. Specified as: a fullword integer; m >= 0.

On Return

c

is the matrix C, containing the coefficients of the spline through the data points (x_j,y_j) for j = 1, n. Returned as: an n by (at least) 4 array, containing numbers of the data type indicated in Table 157.

init

is an indicator that is set to indicate that the coefficients have been initialized. (If you call this subroutine again with the same data, this value should be specified for init.) Returned as: a fullword integer; init = 4.

s

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

Notes

1. In your C program, argument init must be passed by reference.
2. Vectors x, y, and t must have no common elements with matrix C and vector s, and matrix C must have no common elements with vector s; otherwise, results are unpredictable.
3. The elements of vector x must be distinct and must be sorted into ascending order; that is, x₁ < x₂ < ...  < x_n. Otherwise, results are unpredictable. For details on how to do this, see ISORT, SSORT, and DSORT--Sort the Elements of a Sequence.

Function

Interpolation is performed at specified abscissas, t_i for i = 1, m, in vector t, using the cubic spline passing through the data points:

(x_j, y_j)    for j = 1, n

where:

x₁ < x₂ < x₃ < ... < x_n

x_j are elements of vector x.

y_j are elements of vector y.

The value of the cubic spline at each t_i is returned in s_i for i = 1, m. See references [15] and [54]. The coefficients of the spline, c_jk for j = 1, n and k = 1, 4, are returned in matrix 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. The cubic spline values returned in s are computed using the coefficients as follows:

s_i = c_j1 + c_j2 (x_j-t_i) + c_j3 (x_j-t_i)² + c_j4 (x_j-t_i)³    for i = 1, m

where:

j = 1    for t_i <= x₁

j = k    for x₁ < t_i <= x_n, such that x_k-1 < t_i <= x_k

j = n    for x_n < t_i

The values specified for m and init indicate which combination of functions are performed by this subroutine:

• If m = 0 and init > 3, no computation is performed.
• If m = 0 and init <= 3, only the coefficients are computed, and no interpolation is performed.
• If m <> 0 and init > 3, the coefficients are not computed, and the interpolation is performed.
• If m <> 0 and init <= 3, the coefficients are computed, and the interpolation is performed.

In addition, if n = 0, no computation is performed.

The values specified for n and init determine the type of spline function:

• If n = 1, the constructed spline is a constant function.
• If n = 2 and init = 0, the constructed spline is a line through the points.
• If n = 2 and init = 1, the constructed spline is a cubic function through the points whose derivative at x₁ is c₁₁.
• If n = 2 and init = 2, the constructed spline is a cubic function through the points whose derivative at x_n is c₂₁.
• If n = 2 and init = 3, the constructed spline is a cubic function through the points whose derivative at x₁ is c₁₁ and at x_n is c₂₁.

Error Conditions

Computational Errors

None

Input-Argument Errors

1. n < 0
2. m < 0

Example 1

This example computes the spline coefficients through a set of data points with no derivative value specified. It also evaluates the spline at the abscissas specified in T. On output, INIT and C are updated with new values.

Call Statement and Input

             X   Y   C   N  INIT  T   S   M
             |   |   |   |   |    |   |   |
CALL SCSINT( X , Y , C , 6 , 0  , T , S , 4 )

Output

        *                                *
        |  0.000  -0.868   0.000  -0.132 |
        |  1.000  -1.264   0.396  -0.132 |
C    =  |  2.000  -0.076  -1.585   0.660 |
        |  1.100   1.267   0.243  -0.609 |
        |  0.000   1.010   0.014   0.076 |
        | -1.000   0.995   0.000   0.005 |
        *                                *

INIT     =  4
S        =  (-2.792, 1.649, 1.100, -2.000)

Example 2

This example computes the spline coefficients through a set of data points with a derivative value specified at the right endpoint. It also evaluates the spline at the abscissas specified in T. On output, INIT and C are updated with new values.

Call Statement and Input

             X   Y   C   N  INIT  T   S   M
             |   |   |   |   |    |   |   |
CALL SCSINT( X , Y , C , 6 , 2  , T , S , 4 )
 
X        =  (1.000, 2.000, 3.000, 4.000, 5.000, 6.000)
Y        =  (0.000, 1.000, 2.000, 1.100, 0.000, -1.000)

        *              *
        |  .   .  .  . |
        | 0.1  .  .  . |
C    =  |  .   .  .  . |
        |  .   .  .  . |
        |  .   .  .  . |
        |  .   .  .  . |
        *              *
 
T        =  (-1.000, 2.500, 4.000, 7.000)

Output

        *                                *
        |  0.000  -0.865   0.000  -0.135 |
        |  1.000  -1.270   0.405  -0.135 |
C    =  |  2.000  -0.054  -1.621   0.675 |
        |  1.100   1.188   0.379  -0.667 |
        |  0.000   1.303  -0.494   0.291 |
        | -1.000   0.100   1.897  -0.797 |
        *                                *
 
INIT     =  4
S        =  (-2.810, 1.652, 1.100, 1.794)

Example 3

This example computes the spline coefficients through a set of data points with a derivative value specified at both endpoints. It does not evaluate the spline at any points. On output, INIT and C are updated with new values. Because arrays are not needed for arguments t and s, the value 0 is specified in their place.

Call Statement and Input

             X   Y   C   N  INIT  T   S   M
             |   |   |   |   |    |   |   |
CALL SCSINT( X , Y , C , 6 , 3  , 0 , 0 , 0 )
 
X        =  (1.000, 2.000, 3.000, 4.000, 5.000, 6.000)
Y        =  (0.000, 1.000, 2.000, 1.100, 0.000, -1.000)

        *               *
        | -1.0  .  .  . |
        |  0.1  .  .  . |
C    =  |   .   .  .  . |
        |   .   .  .  . |
        |   .   .  .  . |
        |   .   .  .  . |
        *               *

Output

        *                                *
        |  0.000   1.000   3.230   1.230 |
        |  1.000  -1.770  -0.460   1.230 |
C    =  |  2.000   0.079  -1.389   0.310 |
        |  1.100   1.152   0.316  -0.568 |
        |  0.000   1.312  -0.476   0.264 |
        | -1.000  -0.100   1.888  -0.788 |
        *                                *
 
INIT     =  4

Example 4

This example evaluates the spline at a set of points, using the coefficients obtained in Example 3.

Call Statement and Input

             X   Y   C   N  INIT  T   S   M
             |   |   |   |   |    |   |   |
CALL SCSINT( X , Y , C , 6 , 4  , T , S , 4 )

Output