[Quantlib-dev] cubic spline problem

[Ferdinando A. <na...@qu...>]

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