Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I'm in need to implement Monotone Cubic Interpolation for interpolate a sequence of points.
The information I have about the points are x,y and timestamp.

I'm much more an IT guy rather than a mathematical person, so I'm looking for an example of implementation. What I need to do with the resulting functions is store them for future analysis.

My favorite syntax are PHP, Python, Java, Delphi, VB or a generic algorithmic language of your choice.

share|improve this question
    
The wikipedia page you've linked to has links to C++ code at the bottom. –  lhf Jun 14 '11 at 4:06

2 Answers 2

up vote 5 down vote accepted

Not terribly hard to do; as a matter of fact, even if you just limit yourself to Hermite interpolation, there are a number of methods. I shall discuss the three which I have most experience with.

Recall that given points $(x_i,y_i),\quad i=1\dots n$, and assuming no two $x_i$ are the same, one can fit a piecewise cubic Hermite interpolant to the data. (I gave the form of the Hermite cubic in this previous answer.)

To use the notation of that answer, you already have $x_i$ and $y_i$ and require an estimate of the $y_i^\prime$ from the given data. There are at least three schemes for doing this: Fritsch-Carlson, Steffen, and Stineman.

(In the succeeding formulae, I use the notation $h_i=x_{i+1}-x_i$ and $d_i=\frac{y_{i+1}-y_i}{h_i}$.)

The method of Fritsch-Carlson computes a weighted harmonic mean of slopes:

$$y_i^\prime = \begin{cases}3(h_{i-1}+h_i)\left(\frac{2h_i+h_{i-1}}{d_{i-1}}+\frac{h_i+2h_{i-1}}{d_i}\right)^{-1} &\text{ if }\mathrm{sign}(d_{i-1})=\mathrm{sign}(d_i)\\ 0&\text{ if }\mathrm{sign}(d_{i-1})\neq\mathrm{sign}(d_i)\end{cases}$$

the method of Steffen is based on a weighted mean (alternatively, a parabolic fit within the interval):

$$y_i^\prime = (\mathrm{sign}(d_{i-1})+\mathrm{sign}(d_i))\min\left(|d_{i-1}|,|d_i|,\frac12 \frac{h_i d_{i-1}+h_{i-1}d_i}{h_{i-1}+h_i}\right)$$

and the method of Stineman fits to circles:

$$y_i^\prime = \frac{h_{i-1} d_{i-1}h_i^2(1+d_i^2)+h_i d_i h_{i-1}^2(1+d_{i-1}^2)}{h_{i-1} h_i^2(1+d_i^2)+h_i h_{i-1}^2(1+d_{i-1}^2)}$$

The formulae I have given are applicable only to "internal" points; you'll have to consult those papers for the slope formulae for handling the endpoints.


As a demonstration of these three methods, consider these two datasets due to Akima:

$$\begin{array}{|l|l|} \hline x&y\\ \hline 1&10\\2&10\\3&10\\5&10\\6&10\\8&10\\9&10.5\\11&15\\12&50\\14&60\\15&95\\ \hline \end{array}$$

and Fritsch and Carlson:

$$\begin{array}{|l|l|} \hline x&y\\ \hline 7.99&0\\8.09&2.7629\times 10^{-5}\\8.19&4.37498\times 10^{-3}\\8.7&0.169183\\9.2&0.469428\\10&0.94374\\12&0.998636\\15&0.999919\\20&0.999994\\ \hline \end{array}$$

Here are plots of these two datasets:

scatter plots

Here are plots of the cubic spline fits to these two sets:

spline fits

Note the wiggliness that was not present in the original data; this is the price one pays for the second-derivative continuity the cubic spline enjoys.

Here now are plots of interpolants using the three methods mentioned earlier.

This is the Fritsch-Carlson result:

Fritsch-Carlson

This is the Steffen result:

Steffen

This is the Stineman result:

Stineman

(The not-too-good result for the Fritsch-Carlson data set might be due to the use of a cubic Hermite interpolant instead of the rational interpolant Stineman recommended to be used with his derivative prescription.)


As I said, I've had good experience with these three; however, you will have to experiment in your environment on which of these is most suitable to your needs.

share|improve this answer
    
There is a small mistake in the Steffen method: third value in the min() function should be also ABS. $$|\frac{h_i d_{i-1}+h_{i-1}d_i}{h_{i-1}+h_i}|$$ I also can't get the Stineman method to produce smooth results. –  aleksv Nov 16 '13 at 7:45

NB: If I could have made this a comment, I probably would have.

There's a Haskell implementation of natural cubic interpolation, though it doesn't appear to specify monotonicity. It may be worth looking at if you have any background there. If nothing else, it may serve as a starting point. There's also some discussion relative to a 3D rendering engine (which discussion touches on varying interpolation methods, very lightly) here. Getting in touch with some of those involved may lead you to better resources.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.