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.
XKNOT — Array of length NDATA + KORDER containing the knot sequence. (Input)
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 |
|
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. |
3. The spline can be evaluated using BSVAL, and its derivative can be evaluated using BSDER.
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
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
