# Implementation of Monotone Cubic Interpolation

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.

• 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. Jun 14, 2011 at 1:27
• The wikipedia page you've linked to has links to C++ code at the bottom.
– lhf
Jun 14, 2011 at 4:06

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}$$

$$\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:

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

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:

This is the Steffen result:

This is the Stineman result:

(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.

• 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. Nov 16, 2013 at 7:45
• Do this methods produce smooth curves? By sooth I mean with continuous second order derivatives at all points? Apr 3, 2019 at 9:11
• @Confounded, no, only cubic splines have $C^2$ continuity (as mentioned in the answer). In examples like the ones I have shown, you can have $C^2$ or monotonicity, but not both. That's why I showed the cubic spline and also the three monotonic methods for comparison. Apr 3, 2019 at 9:15