cubic spline problem

Posted by Ferdinando Ametrano-3 on Mar 13, 2003; 3:21pm
URL: http://quantlib.414.s1.nabble.com/new-QuantLib-license-tp9975p9976.html

Hi Nicolas

I found out that the algorithm you contributed for cubic spline
interpolation is unable to fit a straight line y(x)=a*x+b when a and b are
not integer. The interpolated function f(x) waves around y(x) and it equals
y(x) at the tabulated x[i],y[i].

As temporary patch I've resumed the old code using the Numerical Recipes
algorithm, and there is a compilation flag (actually a boolean data member)
to switch from the NR algorithm to your code.

Could you take a look at the problem? I don't think it is a round-off
problem, but I might be wrong. I would also appreciate if you could comment
the code a little more, especially when it comes to the definition of the
tridiagonal operator. Specifically I would like to understand which
equation is described in each row of the operator.

I definitely prefer your approach to the NR one, and I hope we could get
the best from both approaches.

Let's try to summarize: we want to interpolate N points using N-1 cubic
polynomials. N-1 cubics are defined by 4(N-1) parameters:
P_i(x)=a_i+b_i*x+c_i*x^2+d_i*x^3.
For each cubic we have 2 equations for the values, P_i(x_i)=y_i and
P_i(x_{i+1})=y_{i+1}, that is 2(N-1) equations, plus 2 equations for the
continuity of the first and second derivatives at the joining knots of the
P_i, that are 2(N-2) equations.
So we have 4N-6 equations and 4N-4 parameters. We need 2 boundary
conditions at X_0 and X_N

NR allows to set at X_0 and X_N the value of the first derivative or to set
the second derivative equal to zero.
If I've got it right your code is smarter in the sense that it guesses the
"best" boundary conditions from the shape of the tabulated (x[i],y[i]), but
it's not clear to me how it does work.

I would like to keep the tridiagonal operator approach (that is more clear
than the cryptic NR code), but we could allow the user to set the first or
second derivative (not necessary equal to zero as in NR) if he likes, and
if he doesn't set them then use your guess algorithm.

Do you agree? Am I missing something? And last but not least would you make
the hacking ;-) ?



------------
ciao -- Nando