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

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

Output

 

KORDER   Maximum difference

 3        0.0096
 4        0.0018
 5        0.0005
 6        0.0004
 7        0.0007
 8        0.0035

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260