I have coded a nice cubic spline interpolator following the basic methods laid out here http://people.math.sfu.ca/~stockie/teaching/macm316/notes/splines.pdf . My program reproduces the example laid out on page 11/15 very well, and I have tested it's ability to interpolate data from many common functions.
I have however, run in to an issue when trying to interpolate data points from the monotonically increasing function $\displaystyle\frac{1}{e^{-0.6x}+1}$.
Image: http://i.stack.imgur.com/JbwoW.png
In the above image the top row is the reproduction of the example from the book that I mentioned. A plot of this function is on the third line. The first plot on the second line is an image of the interpolant function. The image to the right is the first derivative of the interpolant (which is supposed to be continuous using the cubic spline algorithm), and to the right of that, the plot of the second derivative of the interpolant.
Does anyone have a hunch at what is going on here? Thanks!
Best Answer
I assume you are using the "second formulation" from the notes -- you are solving a system of linear equations to get the second derivative vectors $m_i$.
It appears to me that this solving process is not correct in your code. This would explain why your spline has bad (discontinuous) first derivatives. These derivatives are clearly wrong -- regardless of what Matlab or any other package produces, the kind of spline construction you are performing should always produce a curve that has continuous first and second derivatives.
As the graph shows, your spline does have continuous second derivatives. This is to be expected -- no matter what values of $m_i$ you use, the formulae for $a_i$, $b_i$, $c_i$, $d_i$ in step #4 on page 10 will give you a curve that has continuous second derivatives.
There is nothing special about the particular function that you're approximating. I suspect that the other examples you have done are wrong, too, but the errors are too small to see.
What's wrong with your code? My guess is that you have some "off by one" indexing error in the setup of the linear system. But, that's just a guess. If you post your code, we can investigate further.