BSINT

Computes the spline interpolant, returning the B-spline coefficients.

Required Arguments

NDATA — Number of data points. (Input)

XDATA — Array of length NDATA containing the data point abscissas. (Input)

FDATA — Array of length NDATA containing the data point ordinates. (Input)

KORDER — Order of the spline. (Input)

KORDER must be less than or equal to NDATA.

KORDER must be less than or equal to NDATA.

XKNOT — Array of length NDATA + KORDER containing the knot sequence. (Input)

XKNOT must be nondecreasing.

XKNOT must be nondecreasing.

BSCOEF — Array of length NDATA containing the B-spline coefficients. (Output)

FORTRAN 90 Interface

Generic: CALL BSINT (NDATA, XDATA, FDATA, KORDER, XKNOT, BSCOEF)

Specific: The specific interface names are S_BSINT and D_BSINT.

FORTRAN 77 Interface

Single: CALL BSINT (NDATA, XDATA, FDATA, KORDER, XKNOT, BSCOEF)

Double: The double precision name is DBSINT.

Description

Following the notation in de Boor (1978, page 108), let Bj = Bj,k,t denote the j-th B-spline of order k with respect to the knot sequence t. Then, BSINT computes the vector a satisfying

and returns the result in BSCOEF = a. This linear system is banded with at most k − 1 subdiagonals and k − 1 superdiagonals. The matrix

A = (Bj (xi))

is totally positive and is invertible if and only if the diagonal entries are nonzero. The routine BSINT is based on the routine SPLINT by de Boor (1978, page 204).

The routine BSINT produces the coefficients of the B-spline interpolant of order KORDER with knot sequence XKNOT to the data (xi, fi) for i = 1 to NDATA, where x = XDATA and f = FDATA. Let t = XKNOT, k = KORDER, and N = NDATA. First, BSINT sorts the XDATA vector and stores the result in x. The elements of the FDATA vector are permuted appropriately and stored in f, yielding the equivalent data (xi, fi) for i = 1 to N. The following preliminary checks are performed on the data. We verify that

xi < xi+1 | i = 1,…, N - 1 |

ti < ti+1 | i = 1,…, N |

ti ≤ ti+k | i = 1,…, N + k - 1 |

The first test checks to see that the abscissas are distinct. The second and third inequalities verify that a valid knot sequence has been specified.

In order for the interpolation matrix to be nonsingular, we also check tk ≤ x ≤ tN + 1 for i = 1 to N. This first inequality in the last check is necessary since the method used to generate the entries of the interpolation matrix requires that the k possibly nonzero B-splines at xi,

Bj - k +1, …, Bj where j satisfies tj ≤ xi < tj + 1

be well-defined (that is, j − k + 1 ≥ 1).

General conditions are not known for the exact behavior of the error in spline interpolation, however, if t and x are selected properly and the data points arise from the values of a smooth (say Ck) function f, i.e. fi = f(xi), then the error will behave in a predictable fashion. The maximum absolute error satisfies

where

For more information on this problem, see de Boor (1978, Chapter 13) and the references therein. This routine can be used in place of the IMSL routine CSINT by calling BSNAK to obtain the proper knots, then calling BSINT yielding the B-spline coefficients, and finally calling IMSL routine BSCPP to convert to piecewise polynomial form.

Comments

1. Workspace may be explicitly provided, if desired, by use of B2INT/DB2INT. The reference is:

CALL B2INT (NDATA, XDATA, FDATA, KORDER, XKNOT, BSCOEF, WK1, WK2, WK3, IWK)

The additional arguments are as follows:

WK1 — Work array of length (5 * KORDER − 2) * NDATA.

WK2 — Work array of length NDATA.

WK3 — Work array of length NDATA.

IWK — Work array of length NDATA.

2. Informational errors

Type | Code | Description |
---|---|---|

3 | 1 | The interpolation matrix is ill-conditioned. |

4 | 3 | The XDATA values must be distinct. |

4 | 4 | Multiplicity of the knots cannot exceed the order of the spline. |

4 | 5 | The knots must be nondecreasing. |

4 | 15 | The I-th smallest element of the data point array must be greater than the Ith knot and less than the (I + KORDER)-th knot. |

4 | 16 | The largest element of the data point array must be greater than the (NDATA)-th knot and less than or equal to the (NDATA + KORDER)-th knot. |

4 | 17 | The smallest element of the data point array must be greater than or equal to the first knot and less than the (KORDER + 1)st knot. |

Example

In this example, a spline interpolant s, to

is computed. The interpolated values are then compared with the exact function values using the IMSL routine BSVAL.

USE BSINT_INT

USE BSNAK_INT

USE UMACH_INT

USE BSVAL_INT

IMPLICIT NONE

INTEGER KORDER, NDATA, NKNOT

PARAMETER (KORDER=3, NDATA=5, NKNOT=NDATA+KORDER)

!

INTEGER I, NCOEF, NOUT

REAL BSCOEF(NDATA), BT, F, FDATA(NDATA), FLOAT,&

SQRT, X, XDATA(NDATA), XKNOT(NKNOT), XT

INTRINSIC FLOAT, SQRT

! Define function

F(X) = SQRT(X)

! Set up interpolation points

DO 10 I=1, NDATA

XDATA(I) = FLOAT(I-1)/FLOAT(NDATA-1)

FDATA(I) = F(XDATA(I))

10 CONTINUE

! Generate knot sequence

CALL BSNAK (NDATA, XDATA, KORDER, XKNOT)

! Interpolate

CALL BSINT (NDATA, XDATA, FDATA, KORDER, XKNOT, BSCOEF)

! Get output unit number

CALL UMACH (2, NOUT)

! Write heading

WRITE (NOUT,99999)

! Print on a finer grid

NCOEF = NDATA

XT = XDATA(1)

! Evaluate spline

BT = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOEF)

WRITE (NOUT,99998) XT, BT, F(XT) - BT

DO 20 I=2, NDATA

XT = (XDATA(I-1)+XDATA(I))/2.0

! Evaluate spline

BT = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOEF)

WRITE (NOUT,99998) XT, BT, F(XT) - BT

XT = XDATA(I)

! Evaluate spline

BT = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOEF)

WRITE (NOUT,99998) XT, BT, F(XT) - BT

20 CONTINUE

99998 FORMAT (' ', F6.4, 15X, F8.4, 12X, F11.6)

99999 FORMAT (/, 6X, 'X', 19X, 'S(X)', 18X, 'Error', /)

END

Output

X S(X) Error

0.0000 0.0000 0.000000

0.1250 0.2918 0.061781

0.2500 0.5000 0.000000

0.3750 0.6247 -0.012311

0.5000 0.7071 0.000000

0.6250 0.7886 0.002013

0.7500 0.8660 0.000000

0.8750 0.9365 -0.001092

1.0000 1.0000 0.000000