BSOPK
Computes the “optimal” spline knot sequence.
Required Arguments
NDATA — Number of data points. (Input)
XDATA — Array of length NDATA containing the location of the data points. (Input)
KORDER — Order of the spline. (Input)
XKNOT — Array of length NDATA + KORDER containing the knot sequence. (Output)
FORTRAN 90 Interface
Generic: CALL BSOPK (NDATA, XDATA, KORDER, XKNOT)
Specific: The specific interface names are S_BSOPK and D_BSOPK.
FORTRAN 77 Interface
Single: CALL BSOPK (NDATA, XDATA, KORDER, XKNOT)
Double: The double precision name is DBSOPK.
Description
Given the abscissas x = XDATA for an interpolation problem and the order of the spline interpolant k = KORDER, BSOPK returns the knot sequence t = XKNOT that minimizes the constant in the error estimate
∥f − s ∥ ≤ c ∥ f(k) ∥
In the above formula, f is any function in Ck 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, page 204), which in turn is based on a theorem of Micchelli, Rivlin and Winograd (1976).
Comments
1. Workspace may be explicitly provided, if desired, by use of B2OPK/DB2OPK. The reference is:
CALL B2OPK (NDATA, XDATA, KORDER, XKNOT, MAXIT, WK, IWK)
The additional arguments are as follows:
MAXIT — Maximum number of iterations of Newton’s Method. (Input) A suggested value is 10.
WK — Work array of length (NDATA − KORDER) * (3 * KORDER − 2) + 6 * NDATA + 2 * KORDER + 5.
IWK — Work array of length NDATA.
2. Informational errors
|
Type |
Code |
Description |
|
3 |
6 |
Newton’s method iteration did not converge. |
|
4 |
3 |
The XDATA values must be distinct. |
|
4 |
4 |
Interpolation matrix is singular. The XDATA values may be too close together. |
3. The default value for MAXIT is 10, this can be overridden by calling B2OPK/DB2OPK directly with a larger value.
Example
In this example, we compute (for k = 3, …, 8) six spline interpolants sk to F(x) = sin(10x3) on the interval [0, 1]. The routine BSOPK is used to generate the knot sequences for sk and then BSINT is called to obtain the interpolant. We evaluate the absolute error
∣sk − F∣
at 100 equally spaced points and print the maximum error for each k.
USE BSOPK_INT
USE BSINT_INT
USE UMACH_INT
USE BSVAL_INT
IMPLICIT NONE
INTEGER KMAX, KMIN, NDATA
PARAMETER (KMAX=8, KMIN=3, NDATA=20)
!
INTEGER I, K, KORDER, NOUT
REAL ABS, AMAX1, BSCOEF(NDATA), DIF, DIFMAX, F,&
FDATA(NDATA), FLOAT, FT, SIN, ST, T, X, XDATA(NDATA),&
XKNOT(KMAX+NDATA), XT
INTRINSIC ABS, AMAX1, FLOAT, SIN
! Define function and tau function
F(X) = SIN(10.0*X*X*X)
T(X) = 1.0 - X*X
! Set up data
DO 10 I=1, NDATA
XT = FLOAT(I-1)/FLOAT(NDATA-1)
XDATA(I) = T(XT)
FDATA(I) = F(XDATA(I))
10 CONTINUE
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99999)
! Loop over different orders
DO 30 K=KMIN, KMAX
KORDER = K
! Generate knots
CALL BSOPK (NDATA, XDATA, KORDER, XKNOT)
! Interpolate
CALL BSINT (NDATA, XDATA, FDATA, KORDER, XKNOT, BSCOEF)
DIFMAX = 0.0
DO 20 I=1, 100
XT = FLOAT(I-1)/99.0
! Evaluate spline
ST = BSVAL(XT,KORDER,XKNOT,NDATA,BSCOEF)
FT = F(XT)
DIF = ABS(FT-ST)
! Compute maximum difference
DIFMAX = AMAX1(DIF,DIFMAX)
20 CONTINUE
! Print maximum difference
WRITE (NOUT,99998) KORDER, DIFMAX
30 CONTINUE
!
99998 FORMAT (' ', I3, 5X, F9.4)
99999 FORMAT (' KORDER', 5X, 'Maximum difference', /)
END
KORDER Maximum difference
3 0.0096
4 0.0018
5 0.0005
6 0.0004
7 0.0007
8 0.0035
