splineKnots¶

Computes the knots for a spline interpolant

Synopsis¶

splineKnots (xdata)

Required Arguments¶

float xdata[] (Input): Array with ndata components containing the abscissas of the interpolation problem.

Return Value¶

The knots. If the knots cannot be computed, then None is returned.

Optional Arguments¶

order, int (Input)

The order of the spline subspace for which the knots are desired. This option is used to communicate the order of the spline subspace.

Default: order = 4, i.e., cubic splines

opt

This option produces knots that satisfy an optimality criterion.

optItmax, int (Input)

This option allows the user to set the maximum number of iterations of Newton’s method.

Default: optItmax = 10

Description¶

Given the data points x = xdata, the order of the spline k = order, and the number n = ndata of elements in xdata, the default action of splineKnots returns a knot sequence that is appropriate for interpolation of data on x by splines of order k (the default order is \(k = 4\)). The knot sequence is contained in its first n + k positions. If k is even, and we assume that the entries in the input vector x are increasing, then the resulting knot sequence t is returned as

\[\begin{split}\begin{array}{ll} t_i = x_0 & \text{for } i = 0, \ldots, k-1 \\ t_i = x_{j-k/2-1} & \text{for } i = k, \ldots, n-1 \\ t_i = x_{n-1} & \text{for } i = n, \ldots, n+k-1 \\ \end{array}\end{split}\]

There is some discussion concerning this selection of knots in de Boor (1978, p. 211). If k is odd, then t is returned as

\[\begin{split}\begin{array}{c} t_i = x_0 \text{ for } i = 0, \ldots, k - 1 \\ t_i = \left(x_{i - \frac{k-1}{2} - 1} + x_{i - 1 - \frac{k-2}{2}}\right)/2 \text{ for } i = k, \ldots, n - 1 \\ t_i = x_{n-1} \text{ for } i = n, \ldots, n + k - 1 \end{array}\end{split}\]

It is not necessary to sort the values in xdata.

If the option opt is selected, then the knot sequence returned minimizes the constant c in the error estimate

\[∥f - s∥≤c∥f ^{(k)}∥\]

In the above formula, f is any function in \(C^k\), and s is the spline interpolant to f at the abscissas x with knot sequence t.

The algorithm is based on a routine described in de Boor (1978, p. 204), which in turn is based on a theorem of Micchelli et al. (1976).

Examples¶

Example 1¶

In this example, knots for a cubic spline are generated and printed. Notice that the knots are stacked at the endpoints and that the second and next to last data points are not knots.

from numpy import *
from pyimsl.math.splineKnots import splineKnots
from pyimsl.math.writeMatrix import writeMatrix

ndata = 6
xdata = empty(ndata)

for i in range(0, ndata):
    xdata[i] = i
knots = splineKnots(xdata)
writeMatrix("The knots for the cubic spline are:",
            knots, colNumberZero=True)

Output¶

 
                     The knots for the cubic spline are:
          0            1            2            3            4            5
          0            0            0            0            2            3
 
          6            7            8            9
          5            5            5            5

Example 2¶

This is a continuation of the examples for splineInterp. Recall that in these examples, a cubic spline interpolant to a function f is computed first. The values of this spline are then compared with the exact function values. The second example uses a quadratic (\(k = 3\)) and a quintic (\(k = 6\)) spline interpolant to the data. Now, instead of using the default knots, select the “optimal” knots as described above. Notice that the error is actually worse in this case.

from __future__ import print_function
from numpy import *
from pyimsl.math.splineInterp import splineInterp
from pyimsl.math.splineValue import splineValue
from pyimsl.math.splineKnots import splineKnots

# Define function


def F(x):
    return sin(15.0 * x)


# Set up a grid
ndata = 11
fdata = empty((ndata), dtype=double)
xdata = empty((ndata), dtype=double)
for i in range(0, ndata):
    xdata[i] = float(i) / float(ndata - 1)
    fdata[i] = F(xdata[i])

for order in range(3, 7, 3):
    knots = splineKnots(xdata, order=order, opt=True)
    # Compute spline interpolant
    sp = splineInterp(xdata, fdata, order=order,
                      knots=knots)
    # Print results
    print("\nThe order of the spline is: ", order)
    print("     x         F(x)       Interpolant    Error")
    for j in range(int(ndata / 2), int(3 * ndata / 2)):
        x = float(j) / float(2 * ndata - 2)
        y = splineValue(x, sp)
        print("  %6.3f  %10.3f   %10.3f   %10.4f" %
              (x, F(x), y, abs(F(x) - y)))

Output¶

The order of the spline is:  3
     x         F(x)       Interpolant    Error
250      -0.572       -0.543       0.0290
300      -0.978       -0.978       0.0000
350      -0.859       -0.819       0.0401
400      -0.279       -0.279       0.0000
450       0.450        0.429       0.0210
500       0.938        0.938       0.0000
550       0.923        0.879       0.0433
600       0.412        0.412       0.0000
650      -0.320       -0.305       0.0150
700      -0.880       -0.880       0.0000
750      -0.968       -0.920       0.0478

The order of the spline is:  6
     x         F(x)       Interpolant    Error
250      -0.572       -0.578       0.0061
300      -0.978       -0.978       0.0000
350      -0.859       -0.854       0.0054
400      -0.279       -0.279       0.0000
450       0.450        0.448       0.0019
500       0.938        0.938       0.0000
550       0.923        0.920       0.0022
600       0.412        0.412       0.0000
650      -0.320       -0.317       0.0020
700      -0.880       -0.880       0.0000
750      -0.968       -0.966       0.0023

Warning Errors¶

IMSL_NO_CONV_NEWTON Newton’s method iteration did not converge.

Fatal Errors¶

`IMSL_DUPLICATE_XDATA_VALUES`	The xdata values must be distinct.
`IMSL_ILL_COND_LIN_SYS`	Interpolation matrix is singular. The xdata values may be too close together.