Computes the cubic spline interpolant with the ‘not-a-knot' condition.
XDATA — Array of
length NDATA
containing the data point abscissas. (Input)
The data point
abscissas must be distinct.
FDATA — Array of length NDATA containing the data point ordinates. (Input)
BREAK — Array of length NDATA containing the breakpoints for the piecewise cubic representation. (Output)
CSCOEF — Matrix of size 4 by NDATA containing the local coefficients of the cubic pieces. (Output)
NDATA — Number of
data points. (Input)
NDATA must be at least
2.
Default: NDATA = size (XDATA,1).
Generic: CALL CSINT (XDATA, FDATA, BREAK, CSCOEF [,…])
Specific: The specific interface names are S_CSINT and D_CSINT.
Single: CALL CSINT (NDATA, XDATA, FDATA, BREAK, CSCOEF)
Double: The double precision name is DCSINT.
The routine CSINT
computes a C 2 cubic spline
interpolant to a set of data points (xi, fi) for
i = 1,
…, NDATA
= N. The breakpoints of the spline are the abscissas. Endpoint conditions
are automatically determined by the program. These conditions correspond to the
“not-a-knot” condition (see de Boor 1978), which requires that the third
derivative of the spline be continuous at the second and next-to-last
breakpoint. If N is 2 or 3, then the linear or quadratic interpolating
polynomial is computed, respectively.
If the data points arise from the values of a smooth (say C 4) function f, i.e. fi = f(xi), then the error will behave in a predictable fashion. Let ξ be the breakpoint vector for the above spline interpolant. Then, the maximum absolute error satisfies
where
For more details, see de Boor (1978, pages 55− 56).
1. Workspace may be explicitly provided, if desired, by use of C2INT/DC2INT. The reference is:
CALL C2INT (NDATA, XDATA, FDATA, BREAK, CSCOEF, IWK)
The additional argument is
IWK — Work array of length NDATA.
2. The cubic spline can be evaluated using CSVAL; its derivative can be evaluated using CSDER.
3. Note that column NDATA of CSCOEF is used as workspace.
In this example, a cubic spline interpolant to a function F is computed. The values of this spline are then compared with the exact function values.
USE CSINT_INT
USE UMACH_INT
USE CSVAL_INT
IMPLICIT NONE
! Specifications
INTEGER NDATA
PARAMETER (NDATA=11)
!
INTEGER I, NINTV, NOUT
REAL BREAK(NDATA), CSCOEF(4,NDATA), F,&
FDATA(NDATA), FLOAT, SIN, X, XDATA(NDATA)
INTRINSIC FLOAT, SIN
! Define function
F(X) = SIN(15.0*X)
! Set up a grid
DO 10 I=1, NDATA
XDATA(I) = FLOAT(I-1)/FLOAT(NDATA-1)
FDATA(I) = F(XDATA(I))
10 CONTINUE
! Compute cubic spline interpolant
CALL CSINT (XDATA, FDATA, BREAK, CSCOEF)
! Get output unit number.
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99999)
99999 FORMAT (13X, 'X', 9X, 'Interpolant', 5X, 'Error')
NINTV = NDATA - 1
! Print the interpolant and the error
! on a finer grid
DO 20 I=1, 2*NDATA - 1
X = FLOAT(I-1)/FLOAT(2*NDATA-2)
WRITE (NOUT,'(2F15.3,F15.6)') X, CSVAL(X,BREAK,CSCOEF),&
F(X) - CSVAL(X,BREAK,&
CSCOEF)
20 CONTINUE
END
X
Interpolant
Error
0.000
0.000
0.000000
0.050
0.809
-0.127025
0.100
0.997
0.000000
0.150
0.723
0.055214
0.200
0.141
0.000000
0.250
-0.549
-0.022789
0.300
-0.978
0.000000
0.350
-0.843
-0.016246
0.400
-0.279
0.000000
0.450
0.441
0.009348
0.500
0.938
0.000000
0.550
0.903
0.019947
0.600
0.412
0.000000
0.650
-0.315
-0.004895
0.700
-0.880
0.000000
0.750
-0.938
-0.029541
0.800
-0.537
0.000000
0.850
0.148
0.034693
0.900
0.804
0.000000
0.950
1.086
-0.092559
1.000
0.650 0.000000
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |