|
****** |
i n p u t ******
|
|
|
x(1),...,x(npoint) | | |
data abscissae, a s s u m e d to be strictly
increasing .
|
|
y(1),...,y(npoint) | | |
corresponding data ordinates .
|
|
dy(1),...,dy(npoint) | | |
estimate of uncertainty in data, a s s u m-
|
|
npoint.....number of data points, | | |
a s s u m e d .gt. 1
s.....upper bound on the discrete weighted mean square distance of
|
|
|
the approximation | | |
f from the data .
|
|
|
|
|
|
****** |
w o r k a r r a y s *****
v.....of size (npoint,7)
a.....of size (npoint,4)
|
|
|
***** |
o u t p u t *****
a(.,1).....contains the sequence of smoothed ordinates .
|
|
|
a(i,j) = f^(j-1)(x(i)), j=2,3,4, i=1,...,npoint-1 , | | |
i.e., the
|
|
|
first three derivatives of the smoothing spline | | |
f at the
left end of each of the data intervals .
|
|
|
|
w a r n i n g . . . | | |
a would have to be transposed before it
|
|
|
could be used in | | |
ppvalu .
|
|
|
|
|
|
****** |
m e t h o d ******
|
|
|
The matrices | | |
Q-transp*d and Q-transp*D**2*Q are constructed in
|
|
|
|
s e t u p q | | |
from x and dy , as is the vector qty = Q-transp*y .
|
|
|
Then, for given | | |
p , the vector u is determined in c h o l 1 d as
the solution of the linear system
|
|
|
(6(1-p)Q-transp*D**2*Q + p*Q)u | | |
= qty .
|
|
|
|
From |
u , the smoothing spline f (for this choice of smoothing par-
|
|
ameter |
p ) is obtained in the sense that
|
|
|
f(x(.)) |
= y - 6(1-p)D**2*Q*u and
|
|
|
|
f’’(x(.)) |
= 6*p*u .
|
|
|
The smoothing parameter | | |
p is found (if possible) so that
|
|
|
|
with |
sf(p) = s(f) , where f is the smoothing spline as it depends
|
|
on |
p . if s = 0, then p = 1 . if sf(0) .le. s , then p = 0 .
|
|
Otherwise, the secant method is used to locate an appropriate | | |
p in
|
|
the open interval | | |
(0,1) . However, straightforward application of
the secant method, as done in the original version of this program,
|
|
can be very slow and is influenced by the units in which | | |
x and y
are measured, as C. Reinsch has pointed out. Instead, on recommend-
ation from C. Reinsch, the secant method is applied to the function
g:q |--> 1/sqrt{sfq(q)} - 1/sqrt{s} ,
|
|
with |
sfq(q) := sf(q/(1+q)), since 1/sqrt{sfq} is monotone increasing
|
|
and close to linear for larger | | |
q . One starts at q = 0 with a
Newton step, i.e.,
|
|
|
q_0 = 0, | | |
q_1 = -g(0)/g’(0)
|
|
|
|
with |
g’(0) = -(1/2) sfq(0)^{-3/2} dsfq, where dsfq = -12*u-transp*r*u ,
|
|
and |
u as obtained for p = 0 . Iteration terminates as soon as
abs(sf - s) .le. .01*s .
|
|
|
|
|
