|
|
|
from |
* a practical guide to splines * by c. de boor
constructs cubic spline interpolant to given data
tau(i), gtau(i), i=1,...,ntau.
|
|
if |
gamma .gt. 0., 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.
|
|
parameters
input
|
|
tau |
sequence of data points. must be strictly increasing.
|
|
gtau |
corresponding sequence of function values.
|
|
ntau |
number of data points. must be at least 4 .
|
|
gamma |
indicates whether additional flexibility is desired.
= 0., no additional knots
in (0.,3.), 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, i=2,...,ntau-2. see description of meth-
od below. the interpolant gets rounded with increasing
|
|
|
gamma. a value of | | |
2.5 for gamma is typical.
in (3.,6.), same , except that knots might also be added in
intervals in which an inflection point would be permit-
|
|
ted. |
a value of 5.5 for gamma is typical.
output
|
|
|
|
break, coef, l, k |
| |
give the pp-representation of the interpolant.
specifically, for break(i) .le. x .le. 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) and i=1,...,l .
|
|
|
|
iflag |
= 1, ok
= 2, input was incorrect. a printout specifying the mistake
was made.
workspace
|
|
s |
is required, of size (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 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. see below.
|
|
|
|
------ |
m e t h o d ------
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 = z(i) = addg(i+1)/(addg(i) + addg(i+1))
(= .5, in case the denominator vanishes). with
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 .le. z .le. 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) = 2.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 can be computed once fsecnd(i),i=1,...,ntau, is fixed.
|
|
|
|
f |
is automatically continuous and has a continuous second derivat-
ive (except when z = 0 or 1 for some i). we determine fscnd(.) from
|
|
the requirement that also the first derivative of |
| |
f be contin-
uous. 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.
|
|
