CSDEC
Computes the cubic spline interpolant with specified derivative endpoint conditions.
Required Arguments
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)
ILEFT — Type of end condition at the left endpoint. (Input)
ILEFT Condition
0 “Not-a-knot” condition
1 First derivative specified by DLEFT
2 Second derivative specified by DLEFT
DLEFT — Derivative at left endpoint if ILEFT is equal to 1 or 2. (Input)
If ILEFT = 0, then DLEFT is ignored.
IRIGHT — Type of end condition at the right endpoint. (Input)
IRIGHT Condition
0 “Not-a-knot” condition
1 First derivative specified by DRIGHT
2 Second derivative specified by DRIGHT
DRIGHT — Derivative at right endpoint if IRIGHT is equal to 1 or 2. (Input) If IRIGHT = 0 then DRIGHT is ignored.
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)
Optional Arguments
NDATA — Number of data points. (Input)
Default: NDATA = size (XDATA,1).
FORTRAN 90 Interface
Generic: CALL CSDEC (XDATA, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK, CSCOEF [, …])
Specific: The specific interface names are S_CSDEC and D_CSDEC.
FORTRAN 77 Interface
Single: CALL CSDEC (NDATA, XDATA, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK, CSCOEF)
Double: The double precision name is DCSDEC.
Description
The routine CSDEC computes a C2 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 to be selected by the user. The user may specify not-a-knot, first derivative, or second derivative at each endpoint (see de Boor 1978, Chapter 4).
If the data (including the endpoint conditions) arise from the values of a smooth (say C4) 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, Chapter 4 and 5).
Comments
1. Workspace may be explicitly provided, if desired, by use of C2DEC/DC2DEC. The reference is:
CALL C2DEC (NDATA, XDATA, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, 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.
Examples
Example 1
In Example 1, a cubic spline interpolant to a function f is computed. The value of the derivative at the left endpoint and the value of the second derivative at the right endpoint are specified. The values of this spline are then compared with the exact function values.
USE CSDEC_INT
USE UMACH_INT
USE CSVAL_INT
IMPLICIT NONE
INTEGER ILEFT, IRIGHT, NDATA
PARAMETER (ILEFT=1, IRIGHT=2, NDATA=11)
!
INTEGER I, NINTV, NOUT
REAL BREAK(NDATA), COS, CSCOEF(4,NDATA), DLEFT,&
DRIGHT, F, FDATA(NDATA), FLOAT, SIN, X, XDATA(NDATA)
INTRINSIC COS, FLOAT, SIN
! Define function
F(X) = SIN(15.0*X)
! Initialize DLEFT and DRIGHT
DLEFT = 15.0*COS(15.0*0.0)
DRIGHT = -15.0*15.0*SIN(15.0*1.0)
! 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 CSDEC (XDATA, FDATA, ILEFT, DLEFT, IRIGHT, &
DRIGHT, 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 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
Output
X Interpolant Error
0.000 0.000 0.000000
0.050 0.675 0.006332
0.100 0.997 0.000000
0.150 0.759 0.019485
0.200 0.141 0.000000
0.250 -0.558 -0.013227
0.300 -0.978 0.000000
0.350 -0.840 -0.018765
0.400 -0.279 0.000000
0.450 0.440 0.009859
0.500 0.938 0.000000
0.550 0.902 0.020420
0.600 0.412 0.000000
0.650 -0.312 -0.007301
0.700 -0.880 0.000000
0.750 -0.947 -0.020391
0.800 -0.537 0.000000
0.850 0.182 0.000497
0.900 0.804 0.000000
0.950 0.959 0.035074
1.000 0.650 0.000000
Example 2
In Example 2, we compute the natural cubic spline interpolant to a function f by forcing the second derivative of the interpolant to be zero at both endpoints. As in the previous example, we compare the exact function values with the values of the spline.
USE CSDEC_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER ILEFT, IRIGHT, NDATA, NOUT
PARAMETER (ILEFT=2, IRIGHT=2, NDATA=11)
!
INTEGER I, NINTV
REAL BREAK(NDATA), CSCOEF(4,NDATA), DLEFT, DRIGHT,&
F, FDATA(NDATA), FLOAT, SIN, X, XDATA(NDATA), CSVAL
INTRINSIC FLOAT, SIN
! Initialize DLEFT and DRIGHT
DATA DLEFT/0./, DRIGHT/0./
! 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 CSDEC (XDATA, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT,&
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 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
Output
X Interpolant Error
0.000 0.000 0.000000
0.050 0.667 0.015027
0.100 0.997 0.000000
0.150 0.761 0.017156
0.200 0.141 0.000000
0.250 -0.559 -0.012609
0.300 -0.978 0.000000
0.350 -0.840 -0.018907
0.400 -0.279 0.000000
0.450 0.440 0.009812
0.500 0.938 0.000000
0.550 0.902 0.020753
0.600 0.412 0.000000
0.650 -0.311 -0.008586
0.700 -0.880 0.000000
0.750 -0.952 -0.015585
0.800 -0.537 0.000000
