We have seen how it is possible to smoothly interpolate between a set of points \(\left(x_i, y_i\right)\), with the \(x_i\) known as nodes and the \(y_i\) as values, by specifying the gradients \(g_i\) at the nodes and calculating values between adjacent pairs using the uniquely defined cubic polynomials that match the values and gradients at them^{[1]}.
We have also seen how extrapolating such polynomials beyond the first and last nodes can yield less than satisfactory results, which we fixed by specifying the first and last gradients and then adding new first and last nodes to ensure that the first and last polynomials would represent straight lines.
Now we shall see how cubic spline interpolation can break down rather more dramatically and how we might fix it.
Here you can clearly see that even though the values of \(y_i\) are nondecreasing the interpolated values are not. Now this is not just an in issue of aesthetics since there are many kinds of functions that are, by definition, either nondecreasing or nonincreasing, known as monotonic functions. For example, the cumulative distribution function \(F_X\) of a random variable \(X\) must be nondecreasing or else we would be able to find a pair of points \(x_0\) and \(x_1\) such that
A simple way to ensure that such problems don't occur is to insist that each spline is monotonic so that the interpolation \(f\) will be nondecreasing between adjacent pairs of nodes with nondecreasing values and nonincreasing between those with nonincreasing values or, more specifically, that
This moves from a nonmonotonic curve, plotted in red, through a monotonic curve, plotted in green, and back to a nonmonotonic curve as those gradients increase.
The obvious question is whether or not this gives us sufficient control to set the gradients so that every spline is monotonic and, thankfully, there's a very simple way to show that it does; set them all to zero.
Having done so, if the values at a pair of adjacent nodes are equal then the spline connecting them must be a horizontal line since the splines are uniquely defined by the values and gradients at the nodes and a horizontal line meets that definition. If they are not then we need simply note that a cubic polynomial cannot have a derivative of zero at more than two points unless it is a horizontal line and since we have set it to zero at both of the endpoints it can't be zero anywhere between them and so, by the intermediate value theorem, can't change sign.
In both cases the spline is monotonic and so if we set the gradients to zero we will have an interpolation that is monotonic between the nodes, referred to as being piecewise monotonic, although it probably won't be very natural looking, as demonstrated by program 3.
Fortunately for our aesthetic sensibilities we are not limited to setting the endpoint gradients to zero to ensure a cubic spline's monotonicity, as was demonstrated by program 2. Specifically, given left and right endpoints of \(\left(x_0,y_0\right)\) and \(\left(x_1,y_1\right)\), if we define
Otherwise, that they cannot be negative is trivially true since if they were the sign of their associated gradients would be the opposite of that of \(\Delta y\) and the curve would have to change from increasing/decreasing to decreasing/increasing at some point to get from the one endpoint to the other.
To prove that they are bounded above by three we have three shapes of the spline to consider; linear, quadratic and cubic curves.
If the curve is linear then it is trivially monotonic and the endpoint gradients must equal the finite difference \(\frac{\Delta y}{\Delta x}\) and so both \(s_0\) and \(s_1\) will equal one.
Next, if it's quadratic then it must take the form
The turning point \(x^\ast\) at which a quadratic curve
If it equals one then the turning point will occur at infinity and so we may conclude that, if it's quadratic, the curve will be monotonic for all valid values of \(s_0\) and \(s_1\).
Finally, if the curve is cubic then it must take the form
and its turning points will be
The turning points at the extremes of \(s_0\) and \(s_1\) are
Derivation 4 goes further by showing that the curve must be monotonic if either \(s_0\) or \(s_1\) are equal to zero or three.
The easiest way to demonstrate that the conditions for monotonicity are met between the specified bounds is to draw a diagram showing where they are and are not. Program 4 does so by plotting turning points below zero in red, those above one in green and those that don't exist in yellow for \(s_0\) and \(s_1\) between \(1\frac14\) and \(4\frac14\), with \(t^\) at the top and \(t^+\) beneath, marking the bounds with black lines.
Whilst this clearly illustrates that the bounds are sufficient to ensure monotonicity, it also proves that they are not necessary since the yellow ellipses extend beyond them equally for both \(t^\) and \(t^+\). For example, if \(s_0\) equals one and \(s_1\) equals three and a half, we have
We can therefore iterate over the nodes and change their gradients one at a time whenever they would lead to a nonmonotonic spline; if their gradients would need to have different signs to meet the requirements for the splines to their left and to their right then we just set them to zero and, if not, the bound with the smaller magnitude must fall within the bound with the larger and so we simply choose the smaller.
The
Note that we're keeping track of whether any of the gradients were changed and only create a new
Defined in MonoCubicSplineInterpolate.js, the
Program 5 shows how this removes the turning points from the interpolation whilst preserving a relatively natural looking shape that setting the gradients to zero did not.
Whilst this certainly yields an aesthetically pleasing curve, there's no mathematical justification for any particular upper bound on the ratio of the gradients and the finite differences and the choice has a visible effect upon its shape, as demonstrated by program 6.
Fritsch and Carlson recommended constraining the ratios at the endpoints of a spline to a quarter circle of radius three so that
It was because of this that I chose to constrain them to a square that is some fraction of the largest in which the splines are monotonic, although I decided on the default of \(0.89\) because it resulted in a square with an area approximately equal to the quarter circle that they had recommended.
This is indicative of the problem with cubic spline interpolation that I mentioned last time as the reason why I chose not to generalise our implementation to other arithmetic types, like vectors and matrices; the final shape of the interpolation depends upon choices that have nothing whatsoever to do with the points that we are interpolating between. This is in stark contract to linear interpolation which represents, in a mathematically justifiable sense, the best guess at the function that passes through the points.
[2] Akima, H., A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures, Journal of the ACM, Vol. 17, No. 4, 1970.
[3] Fritsch, F. & Carlson, R., Monotone Piecewise Cubic Interpolation, SIAM Journal on Numerical Analysis, Vol. 17, No. 2, 1980.
We have also seen how extrapolating such polynomials beyond the first and last nodes can yield less than satisfactory results, which we fixed by specifying the first and last gradients and then adding new first and last nodes to ensure that the first and last polynomials would represent straight lines.
Now we shall see how cubic spline interpolation can break down rather more dramatically and how we might fix it.
The Problem Of Monotonicity
Program 1 illustrates the problem with the set of points
\[
\left\{(0,10),(2,10),(3,10),(5,10),(6,10),(8,10),(9,10.5),(11,15),(12,50),(14,60),(15,85)\right\}
\]
chosen by Akima^{[2]}.
Program 1: Interpolating Akima's Points



Here you can clearly see that even though the values of \(y_i\) are nondecreasing the interpolated values are not. Now this is not just an in issue of aesthetics since there are many kinds of functions that are, by definition, either nondecreasing or nonincreasing, known as monotonic functions. For example, the cumulative distribution function \(F_X\) of a random variable \(X\) must be nondecreasing or else we would be able to find a pair of points \(x_0\) and \(x_1\) such that
\[
\begin{align*}
x_0 &< x_1\\
F_X\left(x_0\right) &> F_X\left(x_1\right)
\end{align*}
\]
and consequently
\[
\Pr\left(x_0 < X \leqslant x_1\right) = F_X\left(x_1\right)  F_X\left(x_0\right) < 0
\]
which is impossible.A simple way to ensure that such problems don't occur is to insist that each spline is monotonic so that the interpolation \(f\) will be nondecreasing between adjacent pairs of nodes with nondecreasing values and nonincreasing between those with nonincreasing values or, more specifically, that
\[
x_i \leqslant a < b \leqslant x_{i+1} \rightarrow
\begin{cases}
f(a) < f(b) &\quad\quad y_i < y_{i+1}\\
f(a) = f(b) &\quad\quad y_i = y_{i+1}\\
f(a) > f(b) &\quad\quad y_i > y_{i+1}
\end{cases}
\]
where \(\rightarrow\) means implies and the left curly brace is followed by expressions that are conditional upon those to their right.
Tweaking The Gradients
We can exert a lot of control over the shape of a given spline by adjusting the gradients at its endpoints, as demonstrated by program 2.
Program 2: The Effect Of The Endpoint Gradients



This moves from a nonmonotonic curve, plotted in red, through a monotonic curve, plotted in green, and back to a nonmonotonic curve as those gradients increase.
The obvious question is whether or not this gives us sufficient control to set the gradients so that every spline is monotonic and, thankfully, there's a very simple way to show that it does; set them all to zero.
Having done so, if the values at a pair of adjacent nodes are equal then the spline connecting them must be a horizontal line since the splines are uniquely defined by the values and gradients at the nodes and a horizontal line meets that definition. If they are not then we need simply note that a cubic polynomial cannot have a derivative of zero at more than two points unless it is a horizontal line and since we have set it to zero at both of the endpoints it can't be zero anywhere between them and so, by the intermediate value theorem, can't change sign.
In both cases the spline is monotonic and so if we set the gradients to zero we will have an interpolation that is monotonic between the nodes, referred to as being piecewise monotonic, although it probably won't be very natural looking, as demonstrated by program 3.
Program 3: Setting The Gradients To Zero



Fortunately for our aesthetic sensibilities we are not limited to setting the endpoint gradients to zero to ensure a cubic spline's monotonicity, as was demonstrated by program 2. Specifically, given left and right endpoints of \(\left(x_0,y_0\right)\) and \(\left(x_1,y_1\right)\), if we define
\[
\begin{align*}
\Delta x &= x_1  x_0 & \Delta y &= y_1  y_0\\
g_0 &= s_0 \times \tfrac{\Delta y}{\Delta x} & g_1 &= s_1 \times \tfrac{\Delta y}{\Delta x}
\end{align*}
\]
then it will be monotonic if both \(s_0\) and \(s_1\) are no less than zero and no greater than three, as shown by Fritsch and Carlson^{[3]}.
Proving Monotonicity
If \(\Delta y\) equals zero then both \(g_0\) and \(g_1\) will equal zero, as we have already seen is required, for any finite \(s_0\) and \(s_1\), so we might as well just set them both to one in such cases.Otherwise, that they cannot be negative is trivially true since if they were the sign of their associated gradients would be the opposite of that of \(\Delta y\) and the curve would have to change from increasing/decreasing to decreasing/increasing at some point to get from the one endpoint to the other.
To prove that they are bounded above by three we have three shapes of the spline to consider; linear, quadratic and cubic curves.
If the curve is linear then it is trivially monotonic and the endpoint gradients must equal the finite difference \(\frac{\Delta y}{\Delta x}\) and so both \(s_0\) and \(s_1\) will equal one.
Next, if it's quadratic then it must take the form
\[
y = at^2  \frac{s_0}{s_01}at + y_0
\]
for some nonzero value of \(a\), where
\[
t = \frac{xx_0}{\Delta x}\\
\]
so that
\[
\begin{align*}
t_0 &= \frac{x_0x_0}{\Delta x} = \frac{0}{\Delta x} = 0\\
t_1 &= \frac{x_1x_0}{\Delta x} = \frac{\Delta x}{\Delta x} = 1
\end{align*}
\]
and it must be the case that
\[
s_0+s_1=2
\]
as proven in derivation 1.
Derivation 1: The Quadratic Case
Firstly, we have
\[
\begin{align*}
y &= at^2 + bt + c&
y_0 &= at_0^2 + bt_0 + c = c&
y_1 &= at_1^2 + bt_1 + c = a+b+c\\
g &= \frac{\mathrm{d}y}{\mathrm{d}t} \times \frac{\mathrm{d}t}{\mathrm{d}x} = \frac{\mathrm{d}y}{\mathrm{d}x} = \frac{2at+b}{\Delta x}&
g_0 &= \frac{2at_0+b}{\Delta x} = \frac{b}{\Delta x}&
g_1 &= \frac{2at_1+b}{\Delta x} = \frac{2a+b}{\Delta x}
\end{align*}
\]
and therefore
\[
\Delta y = y_1  y_0 = (a+b+c)c = a+b
\]
Recalling that
\[
\begin{align*}
g_0 &= s_0 \times \frac{\Delta y}{\Delta x}&
g_1 &= s_1 \times \frac{\Delta y}{\Delta x}
\end{align*}
\]
we also have the identities
\[
\begin{align*}
g_0 \Delta x &= b = s_0 \Delta y = s_0 (a+b)&
g_1 \Delta x &= 2a+b = s_1 \Delta y = s_1 (a+b)
\end{align*}
\]
Rearranging yields
\[
\begin{align*}
s_0a + \left(s_01\right)b &= 0 & b = \frac{s_0}{s_01}a\\
\left(s_12\right)a + \left(s_11\right)b &= 0 & b = \frac{s_12}{s_11}a
\end{align*}
\]
so that
\[
\begin{align*}
y = at^2  \frac{s_0}{s_01}at + y_0
\end{align*}
\]
and
\[
\begin{align*}
\frac{s_0}{s_01} &= \frac{s_12}{s_11}\\
s_0\left(s_11\right) &= \left(s_01\right)\left(s_12\right)\\
s_0s_1s_0 &= s_0s_12s_0s_1+2\\
s_0 + s_1 &= 2
\end{align*}
\]
The turning point \(x^\ast\) at which a quadratic curve
\[
ax^2+bx+c
\]
changes from increasing to decreasing, or vice versa, occurs when its derivative equals zero
\[
2ax^\ast + b = 0
\]
so that
\[
x^\ast = \frac{b}{2a}
\]
which in this case means that
\[
t^\ast = \frac{s_0}{2s_02}
\]
This can only lie between the nodes of the spline if
\[
0 < \frac{s_0}{2s_02} < 1
\]
Now we know that \(s_0\) cannot be negative and if it is less than one then the denominator will be negative and so
\[
\frac{s_0}{2s_02} \leqslant 0
\]
If it is greater than one then the denominator is positive and we have
\[
\begin{align*}
\frac{s_0}{2s_02} &> 0\\
\frac{s_0}{2s_02} &< 1 \rightarrow
s_0 < 2s_02 \rightarrow
s_0 > 2
\end{align*}
\]
which have already shown cannot be the case since neither \(s_0\) nor \(s_1\) can be negative and their sum must equal two.If it equals one then the turning point will occur at infinity and so we may conclude that, if it's quadratic, the curve will be monotonic for all valid values of \(s_0\) and \(s_1\).
Finally, if the curve is cubic then it must take the form
\[
y = at^3  \frac{2s_0+s_13}{s_0+s_12}at^2 +\frac{s_0}{s_0+s_12}at + y_0
\]
as proven in derivation 2
Derivation 2: The Cubic Case
This time we have
\[
\begin{align*}
y &= at^3 + bt^2 + ct + d&
y_0 &= d&
y_1 &= a+b+c+d\\
g &= \frac{3at^2+2bt+c}{\Delta x}&
g_0 &= \frac{c}{\Delta x}&
g_1 &= \frac{3a+2b+c}{\Delta x}
\end{align*}
\]
and so
\[
\Delta y = y_1y_0 = (a+b+c+d)  d = a+b+c
\]
Multiplying the endpoint gradients by \(\Delta x\) yields
\[
\begin{align*}
g_0 \Delta x &= c = s_0 \Delta y = s_0 (a+b+c)&
g_1 \Delta x &= 3a+2b+c = s_1 \Delta y = s_1 (a+b+c)
\end{align*}
\]
and consequently
\[
\begin{align*}
s_0a + s_0b + \left(s_01\right)c &= 0&
c &= \frac{s_0a + s_0b}{s_01}\\
\left(s_13\right)a + \left(s_12\right)b + \left(s_11\right)c &=0&
c &= \frac{\left(s_13\right)a + \left(s_12\right)b}{s_11}\\
\end{align*}
\]
Solving for \(c\) gives us
\[
\begin{align*}
\frac{s_0a + s_0b}{s_01} &= \frac{\left(s_13\right)a + \left(s_12\right)b}{s_11}\\
\left(\frac{s_0}{s_01}  \frac{s_12}{s_11}\right)b &= \left(\frac{s_13}{s_11}  \frac{s_0}{s_01}\right)a\\
\left(s_0\left(s_11\right)  \left(s_01\right)\left(s_12\right)\right)b &= \left(\left(s_01\right)\left(s_13\right)  s_0\left(s_11\right)\right)a\\
\left(\left(s_0s_1s_0\right)  \left(s_0s_12s_0s_1+2\right)\right)b &= \left(\left(s_0s_13s_0s_1+3\right)  \left(s_0s_1s_0\right)\right)a\\
\left(s_0+s_12\right)b &= \left(2s_0s_1+3\right)a\\
b &= \frac{2s_0+s_13}{s_0+s_12}a
\end{align*}
\]
and so the spline must take the form
\[
\begin{align*}
y &= at^3  \frac{2s_0+s_13}{s_0+s_12}at^2 \frac{s_0a  \frac{2s_0+s_13}{s_0+s_12}s_0a}{s_01}t + y_0\\
&= at^3  \frac{2s_0+s_13}{s_0+s_12}at^2 \frac{\left(\left(s_0+s_12\right)  \left(2s_0+s_13\right)\right)s_0}{\left(s_01\right)\left(s_0+s_12\right)}at + y_0\\
&= at^3  \frac{2s_0+s_13}{s_0+s_12}at^2 +\frac{\left(s_01\right)s_0}{\left(s_01\right)\left(s_0+s_12\right)}at + y_0\\
&= at^3  \frac{2s_0+s_13}{s_0+s_12}at^2 +\frac{s_0}{s_0+s_12}at + y_0
\end{align*}
\]
Note that it cannot be the case that
\[
s_0+s_1=2
\]
since this would lead to infinite coefficients and thus contradict our assumption that the curve is cubic.
and its turning points will be
\[
\begin{align*}
t^ &= \frac{\left(2s_0+s_13\right)  \sqrt{\left(s_0+s_13\right)^2s_0s_1}}{3s_0+3s_16}\\
t^+ &= \frac{\left(2s_0+s_13\right) + \sqrt{\left(s_0+s_13\right)^2s_0s_1}}{3s_0+3s_16}
\end{align*}
\]
as proven in derivation 3.
Derivation 3: The Cubic Turning Points
The derivative of the curve is given by
\[
\frac{\mathrm{d}y}{\mathrm{d}t} = 3at^2  \frac{4s_0+2s_16}{s_0+s_12}at +\frac{s_0}{s_0+s_12}a
\]
and its turning points must therefore satisfy
\[
\begin{align*}
3t^{\ast2}  \frac{4s_0+2s_16}{s_0+s_12}t^\ast +\frac{s_0}{s_0+s_12} &= 0\\
\left(3s_0+3s_16\right)t^{\ast2}  \left(4s_0+2s_16\right)t^\ast + s_0 &= 0
\end{align*}
\]
By the formula for solving quadratic equations we have
\[
\begin{align*}
t^\ast &= \frac{\left(4s_0+2s_16\right) \pm \sqrt{\left(4s_0+2s_16\right)^2  4\left(3s_0+3s_16\right)s_0}}{2\left(3s_0+3s_16\right)}\\
&= \frac{\left(2s_0+s_13\right) \pm \sqrt{\left(2s_0+s_13\right)^2  \left(3s_0+3s_16\right)s_0}}{3s_0+3s_16}\\
&= \frac{\left(2s_0+s_13\right) \pm \sqrt{\left(4s_0^2+s_1^2+9+4s_0s_112s_06s_1\right)  \left(3s_0^2+3s_0s_16s_0\right)}}{3s_0+3s_16}\\
&= \frac{\left(2s_0+s_13\right) \pm \sqrt{s_0^2+s_1^2+9+s_0s_16s_06s_1}}{3s_0+3s_16}\\
&= \frac{\left(2s_0+s_13\right) \pm \sqrt{\left(s_0+s_13\right)^2s_0s_1}}{3s_0+3s_16}
\end{align*}
\]
The turning points at the extremes of \(s_0\) and \(s_1\) are
\[
\begin{matrix}
s_0 & s_1 & t^+ & t^\\
0 & 0 & 0 & 1\\
0 & 3 & 0 & 0\\
3 & 0 & 1 & 1\\
3 & 3 & \frac12 & \frac12\\
\end{matrix}
\]
Note that since we know that the curve is cubic, the last case must be a point of inflection where it goes from increasing to flat to instantaneously back to increasing, or decreasing to flat to decreasing. In the other cases the turning points occur at the endpoints of the curve and so it must be monotonic between them in all four cases.Derivation 4 goes further by showing that the curve must be monotonic if either \(s_0\) or \(s_1\) are equal to zero or three.
Derivation 4: Monotonicity At The Edges
Firstly, if \(s_0\) equals zero then
Secondly, if \(s_0\) equals three we have
Thirdly, if \(s_1\) equals zero then
\[
t^\ast = \frac{\left(s_13\right) \pm \sqrt{\left(s_13\right)^2}}{3s_16}
= \frac{\left(s_13\right) \pm \left(s_13\right)}{3s_16}
\]
and so
\[
\begin{align*}
t^ &= 0&
t^+ &= \frac{2s_16}{3s_16} = \frac{62s_1}{63s_1}
\end{align*}
\]
Since \(t^\) is always zero we require that \(t^+\) cannot be greater than zero or less than one if \(s_1\) lies between zero and three which we can verify with the cases
\[
\begin{align*}
s_0 = 0 &\rightarrow t^+ = 1\\
s_0 < 2 &\rightarrow \left(62s_1>0 \wedge 63s_1>0 \wedge 62s_1 > 63s_1\right) \rightarrow t^+ > 1\\
s_0 = 2 &\rightarrow t^+ = \infty\\
s_0 > 2 &\rightarrow \left(62s_1>0 \wedge 63s_1<0\right) \rightarrow t^+ < 0\\
s_0 = 3 &\rightarrow t^+ = 0
\end{align*}
\]
where \(\wedge\) means and.Secondly, if \(s_0\) equals three we have
\[
t^\ast = \frac{\left(s_1+3\right) \pm \sqrt{s_1^23s_1}}{3s_1+3}
\]
and the term inside the square root will be negative if \(s_1\) is greater than zero or less than three so that the curve won't have any turning points. If it equals zero or three then that term will equal zero and so \(t^\) will equal \(t^+\) and will be a point of inflection.Thirdly, if \(s_1\) equals zero then
\[
t^\ast = \frac{\left(2s_03\right) \pm \sqrt{\left(s_03\right)^2}}{3s_06}
= \frac{\left(2s_03\right) \pm \left(s_03\right)}{3s_06}
\]
and
\[
\begin{align*}
t^ &= \frac{s_0}{3s_06} & t^+ &= 1
\end{align*}
\]
This time \(t^+\) always equals one and we can confirm that \(t^\) lies between zero and one with
\[
\begin{align*}
s_0 = 0 &\rightarrow t^ = 0\\
s_0 < 2 &\rightarrow t^ <0\\
s_0 = 2 &\rightarrow t^ = +\infty\\
s_0 > 2 &\rightarrow t^ > 1\\
s_0 = 3 &\rightarrow t^ = 1
\end{align*}
\]
Finally, \(s_1\) equal to three yields
\[
\begin{align*}
t^\ast &= \frac{2s_0 \pm \sqrt{s_0^23s_0}}{3s_0+3}
\end{align*}
\]
and monotonicity follows by the same argument as made for \(s_0\) equal to three.
The easiest way to demonstrate that the conditions for monotonicity are met between the specified bounds is to draw a diagram showing where they are and are not. Program 4 does so by plotting turning points below zero in red, those above one in green and those that don't exist in yellow for \(s_0\) and \(s_1\) between \(1\frac14\) and \(4\frac14\), with \(t^\) at the top and \(t^+\) beneath, marking the bounds with black lines.
Program 4: The Graphical Solution



Whilst this clearly illustrates that the bounds are sufficient to ensure monotonicity, it also proves that they are not necessary since the yellow ellipses extend beyond them equally for both \(t^\) and \(t^+\). For example, if \(s_0\) equals one and \(s_1\) equals three and a half, we have
\[
\left(s_0+s_13\right)^2s_0s_1
= \left(1+3\tfrac123\right)^21 \times 3\tfrac12
= \left(1\tfrac12\right)^23\tfrac12
= 2\tfrac14  3\tfrac12
= 1\tfrac14
\]
so that the square roots in the formulae for \(t^\) and \(t^+\) are imaginary and consequently the curve doesn't have any turning points.The Piecewise Monotonic Cubic Spline Implementation
A significant benefit of the stricter bounds is that for any \(s_0\) within them every \(s_1\) within them must yield a monotonic curve and viceversa.We can therefore iterate over the nodes and change their gradients one at a time whenever they would lead to a nonmonotonic spline; if their gradients would need to have different signs to meet the requirements for the splines to their left and to their right then we just set them to zero and, if not, the bound with the smaller magnitude must fall within the bound with the larger and so we simply choose the smaller.
The
toMono
function given in listing 1 updates the gradients of the nodes of an ak.cubicSplineInterpolate
interpolator f
so that their ratios to the finite differences between them lie between zero and ub
, ensuring that they fall within the bounds required by the splines to their left and right by iterating over the splines and modifying the gradients at their endpoints so that each will meet the conditions of both.Listing 1: Updating The Gradients
function toMono(ub, f) { var nodes = f.nodes(); var changed = false; var n = nodes.length; var x0 = nodes[0].x; var y0 = nodes[0].y; var g0 = nodes[0].g; var i, x1, y1, g1, g, ubg; for(i=1;i<n;++i) { x1 = nodes[i].x; y1 = nodes[i].y; g1 = nodes[i].g; g = (y1y0)/(x1x0); ubg = ub*g; if(g0!==0) { if(g0*g<=0) {changed=true; nodes[i1].g = 0;} else if(Math.abs(g0)>Math.abs(ubg)) {changed=true; nodes[i1].g = ubg;} } if(g1!==0) { if(g1*g<=0) {changed=true; nodes[i].g = g1 = 0;} else if(Math.abs(g1)>Math.abs(ubg)) {changed=true; nodes[i].g = g1 = ubg;} } x0 = x1; y0 = y1; g0 = g1; } return changed ? ak.cubicSplineInterpolate(nodes) : f; }
Note that we're keeping track of whether any of the gradients were changed and only create a new
ak.cubicSplineInterpolate
interpolator if they were, returning the original one if not.Defined in MonoCubicSplineInterpolate.js, the
ak.monoCubicSplineInterpolate
function given in listing 2 creates a monotonic cubic spline interpolator with an optional first argument, defaulting to \(0.89\) and which we determine has been specified by checking whether or not it's a number, defining how well it should preserve the gradients that an ak.cubicSplineInterpolate
interpolator has after it's been constructed with the remaining arguments. Specifically, the upper bound of the ratio of the gradients to the finite differences passed to toMono
as ub
will be zero if it is zero and three if it's equal to one.Listing 2: ak.monoCubicSplineInterpolate
ak.monoCubicSplineInterpolate = function() { var args = arguments; var arg0 = args[0]; var p = 0.89; var f; if(ak.nativeType(arg0)===ak.NUMBER_T) { p = Number(arg0); if(!(p>=0 && p<=1)) { throw new Error('invalid gradient preservation in ak.monoCubicSplineInterpolate'); } args = Array.prototype.slice.call(args, 1); } try {f = ak.cubicSplineInterpolate.apply(null, args);} catch(e) {throw new Error(e.message.replace('ak.cubic', 'ak.monoCubic'));} return toMono(3*p, f); };
Program 5 shows how this removes the turning points from the interpolation whilst preserving a relatively natural looking shape that setting the gradients to zero did not.
Program 5: Monotonically Interpolating Akima's Points



Whilst this certainly yields an aesthetically pleasing curve, there's no mathematical justification for any particular upper bound on the ratio of the gradients and the finite differences and the choice has a visible effect upon its shape, as demonstrated by program 6.
Program 6: Changing The Upper Bound



Fritsch and Carlson recommended constraining the ratios at the endpoints of a spline to a quarter circle of radius three so that
\[
\begin{align*}
s_0 &\geqslant 0\\
s_1 &\geqslant 0\\
s_0^2+s_1^2 &\leqslant 9
\end{align*}
\]
Unfortunately this introduces a dependency between them which means that the final gradients of the cubic splines will depend upon the order in which we adjust them to meet this condition. That is unless we add a further condition, like choosing the largest ratios that do, which will make them more difficult to calculate.It was because of this that I chose to constrain them to a square that is some fraction of the largest in which the splines are monotonic, although I decided on the default of \(0.89\) because it resulted in a square with an area approximately equal to the quarter circle that they had recommended.
This is indicative of the problem with cubic spline interpolation that I mentioned last time as the reason why I chose not to generalise our implementation to other arithmetic types, like vectors and matrices; the final shape of the interpolation depends upon choices that have nothing whatsoever to do with the points that we are interpolating between. This is in stark contract to linear interpolation which represents, in a mathematically justifiable sense, the best guess at the function that passes through the points.
References
[1] Cubic Line Division, www.thusspakeak.com, 2018.[2] Akima, H., A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures, Journal of the ACM, Vol. 17, No. 4, 1970.
[3] Fritsch, F. & Carlson, R., Monotone Piecewise Cubic Interpolation, SIAM Journal on Numerical Analysis, Vol. 17, No. 2, 1980.
Leave a comment