tautsp Subroutine

If 0 < GAMMA, additional knots are introduced where needed to
make the interpolant more flexible locally.  This avoids extraneous
inflection points typical of cubic spline interpolation at knots to
rapidly changing data.

On the I-th interval, (TAU(I), TAU(I+1)), the interpolant is of the
form:

(*)  F(U(X)) = A + B * U + C * H(U,Z) + D * H(1-U,1-Z),

with

  U = U(X) = ( X - TAU(I) ) / DTAU(I).

Here,

  Z(I) = ADDG(I+1) / ( ADDG(I) + ADDG(I+1) )

but if the denominator vanishes, we set Z(I) = 0.5

Also, we have

  ADDG(J) = abs ( DDG(J) ), 
  DDG(J) = DG(J+1) - DG(J),
  DG(J) = DIVDIF(J) = ( GTAU(J+1) - GTAU(J) ) / DTAU(J)

and

  H(U,Z) = ALPHA * U**3 
         + ( 1 - ALPHA ) * ( max ( ( ( U - ZETA ) / ( 1 - ZETA ) ), 0 )**3

with

  ALPHA(Z) = ( 1 - GAMMA / 3 ) / ZETA
  ZETA(Z) = 1 - GAMMA * min ( ( 1 - Z ), 1/3 )

Thus, for 1/3 <= Z <= 2/3, F is just a cubic polynomial on
the interval I.  Otherwise, it has one additional knot, at

  TAU(I) + ZETA * DTAU(I).

As Z approaches 1, H(.,Z) has an increasingly sharp bend near 1,
thus allowing F to turn rapidly near the additional knot.

In terms of F(J) = GTAU(J) and FSECND(J) = second derivative of F 
at TAU(J), the coefficients for (*) are given as:

  A = F(I) - D
  B = ( F(I+1) - F(I) ) - ( C - D )
  C = FSECND(I+1) * DTAU(I)**2 / HSECND(1,Z)
  D = FSECND(I) * DTAU(I)**2 / HSECND(1,1-Z)

Hence these can be computed once FSECND(1:NTAU) is fixed.

F is automatically continuous and has a continuous second derivative
except when Z=0 or 1 for some I.  We determine FSECND from
the requirement that the first derivative of F be continuous.

In addition, we require that the third derivative be continuous
across TAU(2) and across TAU(NTAU-1).  This leads to a strictly
diagonally dominant tridiagonal linear system for the FSECND(I)
which we solve by Gauss elimination without pivoting.

14 February 2007

Carl DeBoor

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

Input, real ( kind = 8 ) TAU(NTAU), the sequence of data points.  
TAU must be strictly increasing.

Input, real ( kind = 8 ) GTAU(NTAU), the corresponding sequence of
function values.

Input, integer ( kind = 4 ) NTAU, the number of data points.  
NTAU must be at least 4.

Input, real ( kind = 8 ) GAMMA, indicates whether additional flexibility
is desired.
GAMMA = 0.0, no additional knots;
GAMMA in (0.0,3.0), under certain conditions on the given data at
points I-1, I, I+1, and I+2, a knot is added in the I-th interval, 
for I = 2,...,NTAU-2.  See description of method.  The interpolant 
gets rounded with increasing gamma.  A value of 2.5 for GAMMA is typical.
GAMMA in (3.0,6.0), same, except that knots might also be added in
intervals in which an inflection point would be permitted.  A value 
of 5.5 for GAMMA is typical.

Output, real ( kind = 8 ) BREAK(L), real ( kind = 8 ) COEF(K,L), 
integer ( kind = 4 ) L, integer K, give the piecewise polynomial 
representation of the interpolant.  Specifically, 
for BREAK(i) <= X <= BREAK(I+1), the interpolant has the form:
  F(X) = COEF(1,I) +  DX    * ( 
         COEF(2,I) + (DX/2) * (
         COEF(3,I) + (DX/3) *
         COEF(4,I) ) )
with  DX = X - BREAK(I) for I = 1,..., L.

Output, integer ( kind = 4 ) IFLAG, error flag.
1, no error.
2, input was incorrect.

Output, real ( kind = 8 ) S(NTAU,6).  The individual columns of this
array contain the following quantities mentioned in the write up
and below.
S(.,1) = DTAU = TAU(.+1)-TAU;
S(.,2) = DIAG = diagonal in linear system;
S(.,3) = U = upper diagonal in linear system;
S(.,4) = R = right hand side for linear system (initially)
       = FSECND = solution of linear system, namely the second
         derivatives of interpolant at TAU;
S(.,5) = Z = indicator of additional knots;
S(.,6) = 1/HSECND(1,X) with X = Z or 1-Z.