Programming Tips

Mathematica in Education and Research
Vol. 4 No. 3
Copyright 1995 TELOS/Springer-Verlag

Monotonic Interpolation

``The complete assurance that the procedure will never generate 'wild' points makes it attractive as a general-purpose procedure.''
					                             --Russell Stineman, Boeing

by Stan Wagon, Macalester College

Mathematicians and scientists over the centuries have realized how important it is to be able to represent numeric data by a function that passes through the points. Names such as Newton, Lagrange, and Hermite are associated with popular schemes. However, the requirements of the real-world sometimes impose restrictions that make classical approaches inappropriate. BÎzier curves, for example, were born of such real-world necessity (see my column in volume 4, number 1 of this journal). In this note I will discuss an interpolation method developed in 1980 by Russell Stineman, a Boeing employee. It is not at all well known, having appeared in print only in Creative Computing [Sti]; yet it has some important features that are lacking, as far as I know, in all the classical approaches. I present it here so that teachers of numerical analysis might incorporate it into their curriculum and users of interpolation will have access to a robust, well-behaved, and elegant method.

We first show some pitfalls of Mathematica's method, which is a Hermite interpolation with computations based on divided differences. The following data seems to be nicely behaved in that it is increasing and the slopes seem to be steadily decreasing. Yet the interpolating function has a nasty hiccup and two additional unexpected bumps (Figure 1a).

data = {{-1, -1}, {-0.96, -0.4}, {-0.88, 0.3}, 
        {-0.62, 0.78}, {0.13, 0.91}, {1, 1}};

f = Interpolation[data]

Plot[f[t], {t, -1, 1},
  Epilog -> {PointSize[0.025], Point /@ data}];

The method of cubic splines also suffers from unwanted bumpiness, as can be seen in Figure 1b, which shows a spline generated by CubicSpline (in the NumericalMath`SplineFit` package]. This function is a parametric spline function, but the same bumpy behavior is evident with traditional spline functions of the form y = f(x).

<<NumericalMath`SplineFit`

data1 = {{-0.8,-0.9}, {-0.75,0.3},
   {-0.3,0.4}, {0.1,0.5}, {0.5,0.6}};

f = SplineFit[data1, Cubic]

ParametricPlot[f[t], {t, 0, 4},
  Epilog -> {PointSize[0.025], Point /@ data1}];

Figure 1

(a)

(b)

Figure 1.   Monotonic data can lead to nonmonotonic interpolants. These examples show how innocuous-looking data can yield (a) an InterpolatingFunction with three critical points, and (b) a cubic spline interpolant with unwanted bumps.

Code for highly resolved Fig 1 (a)

data = {{-1, -1}, {-0.96, -0.4}, {-0.88, 0.3}, 
        {-0.62, 0.78}, {0.13, 0.91},  {1, 1}};

f = Interpolation[data];
Plot[f[t], {t, -1, 1},
  Epilog->{PointSize[.025], Point /@ data},
  PlotRange->All, DefaultFont->{"Courier", 9},
  MaxBend->1, Ticks->{Range[-1,1,.5],Range[-1,1,.5] }];

Code for highly resolved Fig 1 (b)

<<NumericalMath`SplineFit`

data1 = {{-0.8, -0.9}, {-0.75, 0.3}, 
 {-0.3, 0.4}, {0.1, 0.5}, {0.5, 0.6}};

f = SplineFit[data1, Cubic];

ParametricPlot[f[t], {t, 0, 4},
  Epilog->{PointSize[0.025], Point /@ data1},
  PlotRange->All, Compiled->False,
   MaxBend->1, Ticks->{Range[-.8,.6,.4], Range[-.8,.6,.4] },
  DefaultFont->{"Courier", 9}];

Interpolation does work well on large data sets, however. Users of NDSolve will know that it does a fine job on the numeric data that results from approximate solutions to differential equations. And the resulting InterpolatingFunction objects can be differentiated to check how well the approximate solution fits the differential equation, or plotted to visualize the solution. In my experience, smooth plots always result.

Nevertheless, small data sets do sometimes arise, and that is where a "monotonic method" is called for. We can take the definition of monotonic to be: if ordered x-y data is such that the y-values are monotonic and the slopes of the line segments joining successive points are monotonic, then the interpolating function and its derivative are monotonic; also, if only the successive segment-slopes are monotonic, then the same is true of the derivative. This last clause eliminates unwanted inflection points.

Before describing the Stineman theory, let's look at the results. Figure 2 shows how the Steinman method treats the data of Figure 1(a). The usage is modeled after that of the built-in interpolation.

f = StinemanInterpolation[data]

Plot[f[t], {t,  data[[1, 1]], data[[-1, 1]]},
  Epilog->{PointSize[0.025], Point /@ data}];

Figure 2

Figure 2.   For smoothly increasing data, the Stineman interpolation method produces a smoothly increasing interpolant.

Code for highly resolved Figure 2

data = {{-1, -1}, {-0.96, -0.4}, {-0.88, 0.3}, 
        {-0.62, 0.78}, {0.13, 0.91},  {1, 1}};
        
f = StinemanInterpolation[data];
Plot[f[t], {t,  data[[1, 1]], data[[-1, 1]]},
PlotRange->All, MaxBend->1,
 Ticks->{Range[-1,1,.5],Range[-1,1,.5] },
 Epilog->{PointSize[.025], Point /@ data},
  DefaultFont->{"Courier", 9} ];

The method also works without fuss on larger data sets, such as equispaced points from a sine or logarithm. The error in these cases is a little larger than with Hermite interpolation, but that might be a small price to pay for the robustness that monotonicity guarantees.

We now turn to the details. The basic idea is simple: we first augment each data point (x, y) to a triple (x, y, y') where y' represents the slope at that point. The slope is obtained simply by passing a circle through (x, y) and its two neighboring points on either side and using the slope of the circle at (x, y) to define y'. The endpoints are a special case: if s is the slope of the line from (x1, y1) to (x2, y2), then if s is steeper than y2' we take y1' to be the slope at x1 of a parabola through the two points and having slope y2' at x2; this yields y1' = s + (s— y2'). However, if s is not steeper than y2' we modify this to s + (s— y2')(|s|/(|s| + |s — y2'|)). This ensures that y1' has the same sign as s, which is important for monotonicity. The other endpoint is handled similarly.

The other aspect of the job is to decide how to connect the dots (xi, yi) and (xi+1, yi+1) so as to respect the slopes yi' and yi+1'. Stineman's idea, for a given x between the two x-values, is to start with the linear interpolant ylin. Then let D1 (resp., D2) be the difference between the y-value on the predictor line and ylin; the predictor line refers to the line through (xi, yi) with slope yi' (resp.,  through (xi+1, yi+1) with slope yi+1' ). If D1D2 ‡ 0, then the desired y-value is simply taken to be ylin + (D1 D2)/(D1 + D2). This is simple and has the desired property (at the endpoints the slope agrees with the given slopes). If, however, D1 and D2 have opposite signs, it is a little more complicated; Stineman found that the following expression does the trick:

                                          (D1 D2) (x — x1 + x — x2)
                             ylin  +   -------------------------------   .
                                               (D1 — D2) (x2 — x1)


Figure 3 illustrates the two types of joining functions. It turns out that the inflection point case always yields an interpolant that passes through the midpoint of the line segment containing the two points. We refer the reader to [Sti] for more details of the properties of his joining functions, as well as the algebraic simplifications in the slope computations.

Figure 3

Figure 3.  These diagrams illustrate the two cases that arise in Stineman's choice of functions to join two points with specified slopes.

Code for Figure 3.

pt0 = {0.,.5}; pt1 = {1,1};
s0 = 1.7; s1 = -.7;
f = StinemanInterpolatingFunction[{0,1},
 {{pt0,s0}, {pt1,s1}}];
 
Plot[f[t], {t,0.,1.}, PlotStyle->
{Thickness[.01]},
Epilog-> {Thickness[.0005],

Line[{pt0, pt1}],
Line[{pt0, {pt1[[1]], pt0[[2]]+(pt1[[1]]-pt0[[1]])s0}}],
Line[{pt1, {pt0[[1]], pt1[[2]]+(pt0[[1]]-pt1[[1]])s1}}],
PointSize[.03],
  Point /@ {pt0, pt1}}, Axes->None, MaxBend->1,
  Compiled->False,PlotRange->
  {{-.1, 1.1},{.3, 1.8}}];

pt0 = {0.,.5}; pt1 = {1,1};
s0 = 2; s1 = 1;
f = StinemanInterpolatingFunction[{0,1},
 {{pt0, s0}, {pt1,s1}}];
 
Plot[f[t], {t,0.,1.}, Epilog-> {
Thickness[.0005],
Line[{pt0, pt1}],
Line[{pt0, {pt1[[1]], pt0[[2]]+(pt1[[1]]-pt0[[1]])s0}}],
Line[{pt1, {pt0[[1]], pt1[[2]]+(pt0[[1]]-pt1[[1]])s1}}],
PointSize[.03],
  Point /@ {pt0, pt1}}, Axes->None, MaxBend->1,
  Compiled->False,
  PlotStyle->
{Thickness[.01]},PlotRange->
  {{-.1, 1.1},{.3, 1.5}}];

The accompanying code implements Stineman interpolation. The heart of the routine is compiled for speed. In all other matters, I followed the model of InterpolatingFunction. Note the use of Message[] to return a warning when an argument is outside the interpolation domain. And note how the StinemanInterpolatingFunction object returns a symbolic form unless it is applied to a number, in which case it evaluates. In essence, the Stineman slopes are cached with the function object, and so are not recomputed unnecessarily. A binary search (implemented recursively in home) is used to locate the interval in which a given argument lies. I'm grateful to Rob Knapp of WRI for several helpful suggestions regarding the Stineman code.

I have not included any way of differentiating a StinemanInterpolatingFunction. Implementing that would be a fine exercise! Another approach for users who need that capability would be to wrap a standard Interpolation around an appropriate table of values of the Stineman interpolant; the resulting InterpolatingFunction could then be differentiated.

References

[Sti]   Russell W. Stineman, A consistently well-behaved method of interpolation, Creative  
               Computing, July, 1980, 54—57.

About the Editor

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition and is especially appreciative of the new ways in which Mathematica allows us to look at mathematical objects. He is the author of several books, including Mathematica in Action and, with Ed Packel, Animating Calculus. He also teaches a high-altitude Mathematica course every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalstr.edu

StinemanInterpolationFunction

StinemanInterpolatingFunction::dmval = 
"Input value lies outside domain of the interpolating function.";

Format[StinemanInterpolatingFunction[range_, _]] :=
   SequenceForm["StinemanInterpolatingFunction[", range,", <>]"]

lin = Compile[{x0, y0, x1, y1, x}, (y0(x1-x) + y1(x-x0))/(x1 - x0)];

delta = Compile[{x0, y0, y0p, x, linint}, (y0 + y0p (x - x0) - linint)];

sC = Compile[{x0, x1, x, linint, delta0, delta1},
   linint + delta0 * delta1 * 
    If[delta0 * delta1 > 0, 1/(delta0 + delta1),
     If[delta0 * delta1 == 0, 0, 
        (2x - x0 - x1) / ((delta0 - delta1)(x1 - x0))]]];
    
home[{_, _}, _] := 1;
home[xvals_, x_] := Module[{n = Length[xvals], h},
  h = Quotient[n, 2];
  Which[x < xvals[[h]],   home[Take[xvals, h], x],
        x > xvals[[h+1]], h + home[Take[xvals, h - n], x],
        True, h] ]

endslope[{x1_, y1_}, {x2_, y2_}, s2_] := Module[
  {z = (y2 - y1)/(x2 - x1)}, z + If[!Negative[z (z - s2)], 
      z - s2, 1/(1/(z - s2)  +  Sign[z - s2]/Abs[z])]]

augment[data_] := Module[
 {temp1, temp2, dTrans, yji, ykj, xji, xkj, midslopes},
 dTrans = Transpose[data];
 {xji, yji} = Rest[Drop[#, -1]] - Drop[#, -2] & /@ dTrans;
 {xkj, ykj} = Drop[#, 2] - Rest[Drop[#, -1]]  & /@ dTrans;
 {temp1, temp2} = {xkj^2 + ykj^2, xji^2 + yji^2};
 midslopes = (yji temp1 + ykj temp2) / (xji temp1 + xkj temp2);
 Transpose[{data, 
     Join[{endslope[First[data], data[[2]], First[midslopes]]},
          midslopes,
          {endslope[Last[data], data[[-2]], Last[midslopes]]}]}]]

StinemanInterpolation[data_] := StinemanInterpolatingFunction[
       data[[{1, -1}, 1]], augment[N[data]]]
  
StinemanInterpolatingFunction[{r_, s_}, data_][x_?NumberQ] := 
  Message[StinemanInterpolatingFunction::dmval] /; !(r <= x <= s)
  
StinemanInterpolatingFunction[{r_, s_}, adata_][x_?NumberQ] := 
  Module[{d, linearInt},
   d = Flatten[{
     adata[[home[First /@ First[Transpose[adata]], x] + {0,1}]], x}];
   Apply[sC, 
   Join[d[[{1, 4, 7}]],
        {linearInt = lin @@ d[[{1, 2, 4, 5, 7}]],
         delta @@ Join[d[[{1, 2, 3}]], {x, linearInt}],
         delta @@ Join[d[[{4, 5, 6}]], {x, linearInt}]}  ]]]