cubspl Subroutine

subroutine cubspl(tau, c, n, ibcbeg, ibcend)

************80

! CUBSPL defines an interpolatory cubic spline.

Discussion:

A tridiagonal linear system for the unknown slopes S(I) of
F at TAU(I), I=1,..., N, is generated and then solved by Gauss
elimination, with S(I) ending up in C(2,I), for all I.

Modified:

14 February 2007

Author:

Carl DeBoor

Reference:

Carl DeBoor,
A Practical Guide to Splines,
Springer, 2001,
ISBN: 0387953663,
LC: QA1.A647.v27.

Parameters:

Input, real ( kind = 8 ) TAU(N), the abscissas or X values of
the data points.  The entries of TAU are assumed to be
strictly increasing.

Input, integer ( kind = 4 ) N, the number of data points.  N is
assumed to be at least 2.

Input/output, real ( kind = 8 ) C(4,N).
On input, if IBCBEG or IBCBEG is 1 or 2, then C(2,1)
or C(2,N) should have been set to the desired derivative
values, as described further under IBCBEG and IBCEND.
On output, C contains the polynomial coefficients of
the cubic interpolating spline with interior knots
TAU(2) through TAU(N-1).
In the interval interval (TAU(I), TAU(I+1)), the spline
F is given by
  F(X) = 
    C(1,I) + 
    C(2,I) * H +
    C(3,I) * H**2 / 2 + 
    C(4,I) * H**3 / 6.
where H=X-TAU(I).  The routine PPVALU may be used to
evaluate F or its derivatives from TAU, C, L=N-1,
and K=4.

Input, integer ( kind = 4 ) IBCBEG, IBCEND, boundary condition indicators.
IBCBEG = 0 means no boundary condition at TAU(1) is given.
In this case, the "not-a-knot condition" is used.  That
is, the jump in the third derivative across TAU(2) is
forced to zero.  Thus the first and the second cubic
polynomial pieces are made to coincide.
IBCBEG = 1 means the slope at TAU(1) is to equal the
input value C(2,1).
IBCBEG = 2 means the second derivative at TAU(1) is
to equal C(2,1).
IBCEND = 0, 1, or 2 has analogous meaning concerning the
boundary condition at TAU(N), with the additional
information taken from C(2,N).

Arguments

Type IntentOptional Attributes Name
real(kind=8) :: tau(n)
real(kind=8) :: c(4,n)
integer(kind=4) :: n
integer(kind=4) :: ibcbeg
integer(kind=4) :: ibcend