MATHEMATICA IN ACTION for issue 8-3 of Mathematica in Education and Research

\($DefaultFont\  = \ {"Courier", \ 8}; 

    SetOptions[Plot, \ PlotPoints -> 100]; \)

The Shaped Interpolation Package;
Code For Electronic Version (do not print)

Remove[Explorer`UnusualBases`Harmonic]

BeginPackage["ShapedInterpolation`", "Utilities`FilterOptions`"]; 


Monotonic::usage = "Monotonic is a possible setting to the WeightType option to ShapedInterpolation and InterpolationPlot.";

Convex = "Convex is a possible setting to the WeightType option to ShapedInterpolation and InterpolationPlot.";

AddCircleDerivatives::usage = "AddCircleDerivatives[{{x0,y0}, {x1,y1},...}] uses circles to get derivative data and returns {{x0,{y0,d0}}, {x1,{y1,d1}},...}.";

HermiteRational::usage = "HermiteRational is a possible setting to the WeightType option to ShapedInterpolation and InterpolationPlot.";

Geometric::usage = "Geometric is an option to ShapedInterpolation and InterpolationPlot that specifies that the derivatives should be approximated by a geometric method.";

Arithmetic::usage = "Arithmetic is an option to ShapedInterpolation and InterpolationPlot that specifies that the derivatives should be approximated by an arithmetic method.";

Circles::usage = "Circles is an option to ShapedInterpolation and InterpolationPlot that specifies that the derivatives should be approximated by Stineman's method using the slopes of circles through three consecutive points.";

Harmonic::usage = "Harmonic is an option to ShapedInterpolation and InterpolationPlot that specifies that the derivatives should be approximated by a harmonic method.";

InterpolationPlot::usage = "InterpolationPlot[d] gives a plot of various
interpolating functions through the data set d.";

ShapedInterpolation::usage = "ShapedInterpolation[{{x1, f1}, {x2, f2}, ...
}] constructs an approximate function to interpolate the data. The function
is a shape-preserving rational cubic that attempts to preserve monotonicity
and convexity in the data.";

ShapedInterpolatingFunction::usage = "ShapedInterpolatingFunction[domain, table] represents an approximate function whose values are found by a shape-preserving interpolation method."; 

DerivativeMethod::usage = "DerivativeMethod is an option to InterpolationPlot and ShapedInterpolation that determines the method used to approximate the slopes. Possible settings are Arithmetic, Geometric, and Harmonic.";

WeightType::usage = "WeightType is an option to InterpolationPlot and ShapedInterpolation that determines the weights used in the rational functions. Possible settings are Monotonic, Convex, HermiteRational."; 

PointStyle::usage = "PointStyle is an option to InterpolationPlot that controls the style of the data points.";

ShowPoints::usage = "ShowPoints is an option to InterpolationPlot that controls the display of the data points atop the curve.";

AddDerivativesToDefaultMethod::usage ="AddDerivativesToDefaultMethod is an option to ShapedInterpolation that, provided the method used is the built-in interpolation method (Automatic), adds derivative values obtained by the circle method to the data before calling Interpolation.";

ShapedInterpolatingFunction::nlacng = "The data is convex but the derivative approximations are not interlaced with the slope sequence. This can cause the interpolant to be nonconvex. Try other methods for the DerivativeMethod option."; 

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

Stineman::usage = "Stineman is a possible setting to the Method option of
InterpolationPlot.";

CubicHermite::usage = "CubicHermite is a possible setting to the Method option of InterpolationPlot."; 

Shaped::usage = "Shaped is a possible setting to the Method option of InterpolationPlot.";

StinemanInterpolation::usage = "StinemanInterpolation[{{x1, f1}, {x2, f2}, .. }] constructs am approximate function to interpolate the data. The function is a StinemanInterpolatingFunction object, defined by the Stineman method, which is a monotonic and very well behaved method, though not terribly smooth.";

StinemanInterpolatingFunction::usage = "StinemanInterpolatingFunction[domain, table] represents an approximate function whose values are found by Stineman interpolation.";

NIntegrateWhich::usage = "NIntegrateWhich[e,{x,a,b}] numerically integrates e from a to b where e is a functional expression given by a Which construction.";

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

Begin["`Private`"]; 



(*  The Stineman Method  *)

Format[StinemanInterpolatingFunction[range_, _]] := 
  StringForm["StinemanInterpolatingFunction[``, <>]", range];

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)}, 
   If[z == 0, 0, z + If[Positive[z*(z - s2)], z - s2, 
      1/(1/(z - s2) + Sign[z - s2]/Abs[z])]]]

deriv[data_] := Module[{temp1, temp2, dTrans, yji, ykj, xji, xkj, midslopes}, 
   dTrans = Transpose[data]; {xji, yji} = 
     (Rest[Drop[#1, -1]] - Drop[#1, -2] & ) /@ dTrans; 
    {xkj, ykj} = (Drop[#1, 2] - Rest[Drop[#1, -1]] & ) /@ dTrans; 
    {temp1, temp2} = {xkj^2 + ykj^2, xji^2 + yji^2}; 
    midslopes = (yji*temp1 + ykj*temp2)/(xji*temp1 + xkj*temp2); 
    Join[{endslope[First[data], data[[2]], First[midslopes]]}, midslopes, 
     {endslope[Last[data], data[[-2]], Last[midslopes]]}]]
     
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, 
        (2*x - x0 - x1)/((delta0 - delta1)*(x1 - x0))]]]; 

augment[data_] := Module[
  {temp1, temp2, dTrans, yji, ykj, xji, xkj, midslopes}, 
    dTrans = Transpose[data]; {xji, yji} = 
      (Rest[Drop[#1, -1]] - Drop[#1, -2] & ) /@ dTrans; 
     {xkj, ykj} = (Drop[#1, 2] - Rest[Drop[#1, -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_?MatrixQ] := 
   StinemanInterpolatingFunction[{data[[{1, -1},1]]}, augment[N[data]]]; 

StinemanInterpolatingFunction[{{r_, s_}}, data_][x_?NumberQ] := 
   Message[StinemanInterpolatingFunction::dmval] /;  !r <= x <= s; 

StinemanInterpolatingFunction[{{_, _}}, adata_][x_?NumberQ] := 
   Module[{d, linearInt}, d = 
      Flatten[{adata[[home[First /@ First[Transpose[adata]], x] + {0, 1}]], 
        x}]; 
   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}]}]]
        
sC1 = Compile[{x0, x1, x, linint, delta0, delta1, delta0test, delta1test}, 
    linint + delta0*delta1*If[delta0test*delta1test > 0, 1/(delta0 + delta1), 
       If[delta0test*delta1test == 0, 0, 
        (2*x - x0 - x1)/((delta0 - delta1)*(x1 - x0))]]]; 
        
StinemanInterpolatingFunction[{{_, _}}, adata_][x_Symbol] := 
   Module[{d, linearInt}, 
   ints=Partition[avals=First /@ First /@ adata,2,1];
  
   xmid = Drop[(avals + RotateLeft[avals])/2,-1];
  Off[CompiledFunction::cfr];
ans=Transpose[  {ints, Table[xx = xmid[[i]];
   d = Flatten[{adata[[home[First /@ First[Transpose[adata]], xx] + {0, 1}]], 
        x}]; 
   sC1 @@ 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}],
        delta @@ Join[d[[{1, 2, 3}]], 
             {xx, lin @@ (d[[{1, 2, 4, 5,7}]]/. x->xx)}],  
        delta @@ Join[d[[{4, 5, 6}]], 
              {xx, lin @@ (d[[{1, 2, 4, 5,7}]]/. x->xx)}]}],
        {i, Length[xmid]}]}];
        
  Which @@ Join @@ ({#[[1,1]]<=x<=#[[1,2]],#[[2]]}& /@ ans)]


(* The Main Plotting Function *)

Options[InterpolationPlot] = {Method -> Automatic, ShowPoints -> True, 
 PointStyle -> {PointSize[0.015]}, 
 AddDerivativesToDefaultMethod -> False};

InterpolationPlot[data_, opts___] := 
   Module[{ff, x, mett, ptsQ}, {met, ptsQ, useknotQ, derQ, psty} = 
      {Method, ShowPoints, UseKnots, AddDerivativesToDefaultMethod, 
         PointStyle} /. {opts} /. Options[InterpolationPlot]; 
   
     If[ !ListQ[psty], psty = {psty}];  
     f = Which[(met === Automatic || met === CubicHermite) &&  !derQ, 
       Interpolation[data][x], 
       met === Automatic && derQ, 
       dd = N[deriv[data]]; 
       Interpolation[AddCircleDerivatives[data]][x], 
       met === Shaped, ShapedInterpolation[data, opts][x], 
       met === Stineman, StinemanInterpolation[data][x]]; 
   
     {{x0, x1}, {y0, y1}} = 
      {({Min[#], Max[#]} & ) @ (First /@ data), 
       ({Min[#], Max[#]} & ) @ (Last /@ data)}; 
    
     prd = PlotRange[Plot[Evaluate[f], {x, x0, x1}, 
            Evaluate[FilterOptions[Plot, opts]],
            PlotRange->All,
            DisplayFunction->Identity]];

      y0 = Min[y0, prd[[2,1]]];
      y1 = Max[y1, prd[[2,2]]];
       
      xr = x1-x0; yr = y1-y0;
      x0 -= xr/20; x1 += xr/20; y0 -= yr/15; y1 += yr/15;
           
     If[ptsQ, Show[Block[{$DisplayFunction = Identity}, 
        (Plot[Evaluate[f], {x, #1[[1]], #1[[2]]}, 
            Evaluate[FilterOptions[Plot, opts]]] & ) /@ 
         Partition[First /@ data, 2, 1]], 
       Evaluate[FilterOptions[Graphics, opts]], 
      
       Frame -> True, 
       FrameTicks->{Automatic,Automatic,None,None},
       Axes -> None, 
       PlotRange -> {{x0, x1}, {y0, y1}}, 
       Epilog -> (jj=Join[psty, Point /@ data])], 
      Plot[Evaluate[f], {x, f[[0,1,1,1]], f[[0,1,1,2]]}, 
       Evaluate[FilterOptions[Plot, opts]], Frame -> True, FrameTicks->{Automatic,Automatic,None,None},
       Axes -> None, 
       PlotRange -> {{x0, x1}, {y0, y1}}, 
       Epilog -> If[ptsQ, {PointSize[0.02], Point /@ data}, {}]]]];

ShapedInterpolatingFunction[{r_, s_}, data_][(x_)?NumberQ] := 
   Message[ShapedInterpolatingFunction::dmval] /;  !r <= x <= s; 

Format[ShapedInterpolatingFunction[range_, _]] := 
   StringForm["ShapedInterpolatingFunction[``, <>]", range]; 
   
  diffs[data_, pos_] := 
   Apply[Subtract, Reverse /@ Partition[Transpose[N[data]][[pos]], 2, 1], 
    {1}]; 

Options[ShapedInterpolation] := {DerivativeMethod -> Automatic, 
   WeightType -> Automatic}; 

AddCircleDerivatives[data_] :=  (dd = N[deriv[data]]; 
    Transpose[{First /@ data, 
            Transpose[{Last /@ data, dd}]}]) 

 ShapedInterpolation[data_?MatrixQ, opts___] :=
  Module[{n=Length[data], xvals, yvals, int, α, Δ, dd, rr}, 
     {met, type} = {DerivativeMethod, 
      WeightType} /. {opts} /. 
       Options[ShapedInterpolation]; 
   
     xvals = First /@ N[data]; yvals = Last /@ N[data]; 
     monoQ = monotonicQ[data]; conQ = convexQ[data]; 
     If[type === Automatic && monoQ && !conQ, type = Monotonic]; 
     If[type === Automatic && conQ, type = Convex]; 
     If[type === Automatic && !monoQ && !conQ, type = Convex]; 
     If[met === Automatic && monoQ &&  !conQ, met = Geometric]; 
     If[met === Automatic &&  !monoQ && conQ, met = Arithmetic]; 
     If[met === Automatic && monoQ && conQ, met = Harmonic]; 
     If[met === Automatic &&  !monoQ &&  !conQ, met = Circles]; 
     int[1] = {2, 3}; int[n] = {n - 2, n - 1}; int[i_] := {i - 1, i + 1}; 
     Δ[i_] := Δ[i, i + 1]; 
     Δ[i_, j_] := 
      (yvals[[j]] - yvals[[i]])/(xvals[[j]] - xvals[[i]]); 
     α[i_, j_] := 
      Module[{xk = xvals[[DeleteCases[int[i], j]]]}, 
       Times @@ ((xk - xvals[[i]])/(xk - xvals[[j]]))]; 
  d = Which[
     met === Geometric, 
       Table[Product[If[(bb = Δ[i, int[i][[k]]]) == 0, 0, 
          Sign[bb]*Abs[bb]^α[i, int[i][[k]]]], {k, Length[int[i]]}], 
        {i, n}], 
     met === Arithmetic, 
       Table[Sum[α[i, int[i][[k]]]*Δ[i, int[i][[k]]], 
         {k, Length[int[i]]}], {i, n}], 
     met === Harmonic, 
       Table[Chop[1/
          Sum[α[i, int[i][[k]]]*
            If[(dd = Δ[i, int[i][[k]]]) == 0, 10^100, 1/dd], 
           {k, Length[int[i]]}]], {i, n}],
     met === Circles, N[deriv[data]], 
     met === Automatic, {}]; 
    
     r = Which[met === Automatic || met === CubicHermite, Array[3 &, {n}], 
       type === HermiteRational, Table[3, {n}], 
       type === Monotonic, 
       Chop[Table[If[(denom = Δ[i]) == 0, 3, 
          1 + (d[[i]] + d[[i + 1]])/denom], {i, n - 1}]], 
       type === Convex1, 
       2 + Table[(mm = 
            Max[d[[i + 1]] - Δ[i], 
             Δ[i] - d[[i]]])/
          (Min[d[[i + 1]] - Δ[i], 
             Δ[i] - d[[i]]] /. 0 | 0. -> mm), {i, n - 1}], 
       type === Convex, Table[1 + 
         (d[[i + 1]] - Δ[i])/
          (Δ[i] - d[[i]] /. 
            0 | 0. -> (d[[i + 1]] - Δ[i])/3) + 
         (Δ[i] - d[[i]])/
          (d[[i + 1]] - Δ[i] /. 
            0 | 0. -> (Δ[i] - d[[i]])/3), {i, n - 1}]]; 
     ShapedInterpolatingFunction[{data[[{1, -1},1]]}, 
      {xvals, yvals, diffs[data, 1], d, r}]];

ShapedInterpolatingFunction[{{_, _}}, intdata_][x_?NumberQ] := 
  Module[{xvals, yvals, dx, d, r, i, t}, 
   xvals = intdata[[1]]; 
   yvals = intdata[[2]]; 
   dx = intdata[[3]]; 
   d = intdata[[4]];
   r = intdata[[5]];
   i = home[xvals, x]; 
   t = (x - xvals[[i]])/dx[[i]];
 
   If[Global`traceQ, AppendTo[Global`ans, 
   {r[[i]], (Global`t-xvals[[i]])/dx[[i]], yvals[[i]], yvals[[i + 1]], dx[[i]], d[[i]], 
      d[[i + 1]]}]];
      
       Return[ratCom[r[[i]], t, yvals[[i]], yvals[[i + 1]], dx[[i]], d[[i]], 
      d[[i + 1]]]]; (* old version: 
   {yvals[[i + 1]], 
       r[[i]]*yvals[[i + 1]] - dx[[i]]*d[[i + 1]], 
       r[[i]]*yvals[[i]] + dx[[i]]*d[[i]], yvals[[i]]} . 
      {t^3, t^2*(1 - t), t*(1 - t)^2, (1 - t)^3}/(1 + (r[[i]] - 3)t(1 - t))*)];
      
      
ShapedInterpolatingFunction[{{_, _}}, intdata_][x_Symbol] := 
  Module[{xvals, yvals, dx, d, r, i, t}, 
 
  Off[CompiledFunction::cfr];
   xvals = intdata[[1]]; 
   xmid = Drop[(xvals + RotateLeft[xvals])/2.,-1];
   yvals = intdata[[2]]; 
   dx = intdata[[3]]; 
   d = intdata[[4]];
   r = intdata[[5]];
ans= Transpose[{  Table[xx=xmid[[j]];
   i = home[xvals, xx]; 
  
   t = (x - xvals[[i]])/dx[[i]];
 Simplify@Chop@Expand@  ratCom[r[[i]], t, yvals[[i]], yvals[[i + 1]], dx[[i]], d[[i]], 
      d[[i + 1]]],{j, Length[xmid]}], Partition[xvals, 2, 1]}];
      On[CompiledFunction::cfr];
     Which @@ Flatten[Join[{(#[[1,1]] <= x <= #[[1,2]]), #[[2]]}& /@ (Reverse /@ ans)],
    1]
       ]
      
 ratCom = Compile[{rr, tt, y1, y2, dx1, d1, d2}, 
    (y2*tt^3 + (rr*y2 - dx1*d2)*tt^2*(1 - tt) + 
       (rr*y1 + dx1*d1)*tt*(1 - tt)^2 + y1*(1 - tt)^3)/
     (1 + (rr - 3)*tt*(1 - tt))]; 




(* Tests for Monotonicity and Convexity *)

monotonicQ[yvals_?VectorQ] := 
  Length[DeleteCases[Union[Sign[Apply[Subtract, 
        Reverse /@ Partition[N[yvals], 2, 1], {1}]]], 0]] != 2;

monotonicQ[data_?MatrixQ] := monotonicQ[Last /@ data];

convexQ[data_] := 
  Module[{xvals, yvals}, {xvals, yvals} = Transpose[data]; 
    monotonicQ[Apply[Subtract, 
      Reverse /@ Partition[yvals, 2, 1], {1}]/
      Apply[Subtract, Reverse /@ Partition[xvals, 2, 1], {1}]]];

NIntegrateWhich[s_, {x_, a_, b_}]:= Module[{d},
  Plus @@( 
  (d = #[[2]];
  If[d === 0 || d===0.||Max[a,#[[1,1]]]>=Min[b, #[[1,2]]],0,
		NIntegrate[d, {x, Max[a,#[[1,1]]],Min[b, #[[1,2]]]}]])& /@ 
({{#[[1,1]], #[[1,3]]}, #[[2]]}& /@ Partition[List@@s,2,2]))]

End[]; EndPackage[];

Monotonic Interpolation, Revisited

Getting a smooth curve through some data points is not always enough. Sometimes it is critical to preserve the general shape suggested by the data.

by  Stan Wagon

.  Introduction

In a column four years ago [Wagon 1995] I discussed the Stineman interpolation method, which, given a sequence of monotonically increasing points, produces a monotonic interpolant. I have since learned about other methods developed about the same time as Stineman’s that yield somewhat better results, and I will discuss one of those here. The problem of monotonic interpolation is an important one in computer-aided design, and this explains why many of the workers in the field are in industry or departments of engineering. However, the issue of monotonic interpolation arises in other contexts as well. Two examples: (1) One might have some wind-tunnel data, which can be a very sparse collection of data. The coefficient–angle relationship is known to be monotonic, and one wants the interpolating curve to respect that.  (2) One might be making an approximation to a probability distribution (Ralph Miles of JPL observes that these arise in the context of risk assessment); the distribution is monotonic and one might have only a very few data points.

The technical definition is as follows: A method that produces an interpolating function for data points ((Subscript[x, i]), (Subscript[y, i])) in the plane is called monotonic if (1) whenever the data are monotonic, so is the function, and (2) whenever the slopes of the data are monotonic, so are the slopes of the interpolant. Condition (2) is also called preservation of convexity. The slopes of the data refer to the slopes of the lines connecting neighboring points.

To set the stage, we recall how the built-in interpolation method (by default it uses a cubic Hermite formula) and the Stineman method (a rational function) do on some monotonic data.

Here are some sample data sets from the literature that will serve as examples.

\(data1
       = {{1, 10}, {2, 10}, {3, 10}, {5, 10}, {6, 10}, {8, 10}, {9, 10.5}, {
          11, 15}, {12, 50}, {14, 60}, {15, 95}}; 


    data2 = {{0, 10}, {1, 4}, {4, 2}, {6.5, 1}, {9, 3}, {10, 10}}; 


    data3 = {{7.99, 0}, {8.09, 0.000027629}, {8.19, 0.00437498}, {8.7, 
          0.169183}, {9.2, 0.469428}, {10, 0.94374}, {12, 0.998636}, {15, 
          0.999919}, {20, 0.999994}}; 


    data4 = {{22, 523}, {22.5, 543}, {22.6, 550}, {22.7, 557}, {22.8, 565}, {
          22.9, 575}, {23, 590}, {23.1, 620}, {23.2, 860}, {23.3, 915}, {
          23.4, 944}, {23.5, 958}, {24, 986}}; 


    data5 = {{0, 10.1}, {1.1, 8}, {1.31, 4.7}, {2.5, 4.}, {3.9, 3.48}, {4.4, 
          3.3}, {5.5, 5.8}, {6, 7}, {8, 7.7}, {10.1, 8.6}}; \)

Recall that Interpolation[data] gives an InterpolatingFunction object that represents a function that interpolates the data set.

Interpolation[data1]

We could now use Plot, but I will instead use a specialized routine from the ShapedInterpolation package included in the electronic supplement to this article; that allows one to specify the method of interpolation, and the resulting plot shows the points and the curve. The default method is Mathematica’s standard Interpolation. As can be seen in Figure , the default method introduces a bump: it is not a monotonic method. This sort of behavior is also typical of many spline methods.

InterpolationPlot[data1];

Figure  shows Stineman’s monotonic method [Stineman 1980] on the same data set. The result is indeed monotonic, but the curve is not as well shaped as one might hope: despite the fact that the first derivatives match at the knots, the changes seem a little abrupt.

InterpolationPlot[data1, Method -> Stineman, 
      ShowPoints -> False];

.  Shaplier Interpolation

We will briefly discuss some methods due to Delbourgo and Gregory [Delbourgo and Gregory 1985] that yield shaplier curves than the Stineman method. The main algorithm in this article is a cubic version that I learned from Brenda Kroschel [Anderson 1987]; it is a variation of the quadratic method of Delbourgo and Gregory.

Suppose we are given data ((Subscript[x, i]), (Subscript[y, i])) with (Subscript[x, i]) < (Subscript[x, i + 1]). First define (Subscript[d, i]) to be an approximation to the derivative at (Subscript[x, i]) of the function that we assume underlies the data. Stineman’s method also needs such values; he used circles to generate them. We will present some more sophisticated ways of getting the d-values in a moment. Define auxiliary terms as follows:

(Subscript[h, i])=(Subscript[x, i + 1]) - (Subscript[x, i])      \(TraditionalForm
      \`Δ\_i = \(y\_\(i + 1\) - y\_i\)\/h\_i\)      θ(x) = ((x - (Subscript[x, i]))/(Subscript[h, i]))  (we often use just θ)

Then if we define polynomials P and Q as follows and use appropriate coefficients (Subscript[r, i]), their quotient s(x) will have some nice properties.

\(TraditionalForm
      \`\(P\_i\)(θ) = 
        \(y\_\(i + 1\)\) θ\^3 + 
          \((\(r\_i\) y\_\(i + 1\) - \(h\_i\) d\_\(i + 1\))\) 
            \(\(θ\^2\)(1 - θ)\) + 
          \((\(r\_i\) y\_i + \(h\_i\) d\_i)\) \(θ(1 - θ)\)\^2 + 
          \(\(y\_i\)(1 - θ)\)\^3\)

          \(TraditionalForm
      \`\(Q\_i\)(θ) = 
        θ\^3 + 
          \(\(r\_i\)(
            \(θ\^2\)(1 - θ) + \(θ(1 - θ)\)\^2 + 
                \((1 - θ)\)\^3 = 
              1 + \((r\_i - 3)\) \(θ(1 - θ)\)\)\)

          s(x) = (((Subscript[P, i])(θ))/((Subscript[Q, i])(θ))) where x is between (Subscript[x, i]) and (Subscript[x, i + 1])

The interpolating function s will pass through the data, and have its derivative at (Subscript[x, i]) equal (Subscript[d, i]). And, under most conditions, the monotonicity and/or convexity of the data will be respected. We say “most” because there are issues having to do with how (Subscript[r, i]) and (Subscript[d, i]) are defined. Here are the three main methods of getting (Subscript[d, i]) (implemented in the package as options). The formulas below apply when 1 < i < n; slightly different conditions are used at the two endpoints (extrapolation based on (Subscript[Δ, 2]) and (Subscript[Δ, 3]), or (Subscript[Δ, n - 2]) and (Subscript[Δ, n - 1])).

• Arithmetic:  \(TraditionalForm
      \`d\_i = \(\(h\_i\) Δ\_\(i - 1\) + 
            \(h\_\(i - 1\)\) Δ\_i\)\/\(h\_i + 
            h\_\(i - 1\)\)\).

• Geometric:  \(TraditionalForm
      \`d\_i = \((
            \(Δ\_\(i - 1\)\^h\_i\) 
              Δ\_i\^h\_\(i - 1\))\)\^\(1/
            \((h\_\(i - 1\) + h\_i)\)\)\)

• Harmonic:  (Subscript[d, i])=(((Subscript[h, i]) + 
                  (Subscript[h, i - 1]))/((Subscript[h, i]) (1/(Subscript[Δ, i - 1]))  
                   +   (Subscript[h, i - 1]) (1/(Subscript[Δ, i]))))

And here is how the r-values can be defined. The best choice depends on whether the data are monotonic or convex or both, and that can be detected by the program.

• monotonic, but not convex:  \(TraditionalForm
      \`r\_i = 1 + \(d\_i + d\_\(i + 1\)\)\/Δ\_i\)

• otherwise:  \(TraditionalForm
      \`r\_i = 1 + 
          \(d\_\(i + 1\) - Δ\_i\)\/\(Δ\_i - 
              d\_i\) + 
          \(Δ\_i - d\_i\)\/\(d\_\(i + 1\) - 
              Δ\_i\)\)

• (Subscript[r, i]) = 3 is also a possibility, and yields a standard cubic Hermite rational interpolation based on the data and derivative values.

Here is an example. By using the Shaped option, the cubic method just described will be used. The package detects the monotonicity and convexity in the data and uses the appropriate r-formula. And it also uses this monotonic/convex information to choose one of the three derivative methods according to recommended practice:

     monotonic and nonconvex -> Geometric
     convex and nonmonotonic -> Arithmetic
     monotonic and convex -> Harmonic
     nonmonotonic and nonconvex -> Circles (the Stineman method)

InterpolationPlot[data1, Method -> Shaped, 
      ShowPoints -> False];

For ease of comparison, Figure  compares the new method with the Stineman method, focusing on the region where the methods differ.

\(\(Show[
      Block[{$DisplayFunction = Identity}, {
			
          InterpolationPlot[data1, Method -> Stineman, 
				\ \ \ 
            ShowPoints -> False, 
			\ \ \ 	
            PlotStyle -> {GrayLevel[0.6], Thickness[0.008]}], 
			
          InterpolationPlot[data1, Method -> Shaped, 
	\ \ \ 			
            ShowPoints -> False, 
			\ \ \ 	
            PlotStyle -> {Thickness[0.002]}]}], 
	\ 
      PlotRange -> {{9.5, 14.5}, {6, 80}}]; \)\)

Figure  shows three more examples using some data sets from the literature. In all cases the results using the shaped interpolant are good (though it must be admitted that sometimes the Stineman method looks better on parts of the domain). The reader can try these with Method set to Automatic, which causes the default interpolation to be used; in the first two cases this leads to unwanted bumps and cusps.

Code for Figure  (Do not print)

\(\(\(Show[
      GraphicsArray[
			
        Block[{$DisplayFunction = Identity}, {

            InterpolationPlot[data2, Method -> Shaped, 
		
              PointStyle -> PointSize[0.02], 
					
              FrameTicks -> {Automatic, \ Range[2, 8, 3], None, None}], 

            InterpolationPlot[data3, Method -> Shaped, 
		
              PointStyle -> PointSize[0.02], 
					
              FrameTicks -> {Automatic, \ Range[0, 1, .25], None, None}], 

            InterpolationPlot[data4, Method -> Shaped, 
					
              FrameTicks -> {Automatic, \ Range[500, 900, 200], None, None}, 
              
	PointStyle -> PointSize[0.02]]}]]]; \)
	\)\)

.  Misshapen Data

Some of these ideas are relevant to the problem of interpolating data that is neither monotonic nor convex, but for which there is a rough shape that one wants to preserve. Keep in mind that the literature on this subject is quite large, making use of ideas of splines. But we present a simple method that can be helpful. The ShapedInterpolation package allows one to use the built-in interpolation method, but first adding derivative values to the data. The circle method of getting derivative values seems to be a good one, and that is what is called when the AddDerivativesToDefaultMethod option is used. First observe (Figure ) that the basic interpolation method is not so good.

InterpolationPlot[data5, PlotRange -> {-3, 15}, 
      Method -> Automatic, PointStyle -> PointSize[0.02],  
      FrameTicks -> {Automatic,  Range[0, 15, 3],  None,  None}];

But when some first-derivative values are added to the data (this works because the built-in Interpolation can handle data points of the form {x, {y, (y^′), …}}), the result, shown in Figure , is quite good (and smoother-looking than Stineman’s method).

\(\(InterpolationPlot[data5, Method -> Automatic, 
      AddDerivativesToDefaultMethod -> True, 
	
      PointStyle -> PointSize[0.02], 
	PlotRange -> {\(-3\), 15}, 
      
\ \ FrameTicks 
        -> {Automatic, \ Range[0, 15, 3], \ None, \ None}]; \)\)

The function of Figure  can also be obtained directly by using the package function AddCircleDerivatives as follows.

newdata5 = AddCircleDerivatives[data5]

Interpolation[newdata5]

.  Implementation and an Application

The implementation of this new method is much the same as the Stineman implementation described in my earlier column. A ShapedInterpolatingFunction object is introduced, with a new feature being that it can accept options. This object contains various numbers important to the computation, and the computation then gets carried out when a numerical argument is given. The function object  is created by calling ShapedInterpolation on a data set. Compilation is used as appropriate and the object is formatted the same was as the built-in InterpolatingFunction objects are, via a line as follows:

\(\(\(Format[ShapedInterpolatingFunction[range_, \ _]]\  := \ 
\ \ \ 
      StringForm["ShapedInterpolatingFunction[``, <>]", \ range]; \)\ 
    \)\)

One new feature of the present work is that both ShapedInterpolatingFunctions and StinemanInterpolatingFunctions can be evaluated for symbolic inputs. This was done in the package by feeding the symbol, say x, to the various routines that do the rational arithmetic that produces the function. These are compiled, so a certain compilation error message is turned off and the compiled functions just use standard evaluation to deal with the nonnumeric argument x. The function is returned as a Which statement that deals with all the intervals. One application of this is that we can easily look at the derivatives (and integrals) of the shaped interpolating functions. To illustrate, let us look at the two examples involving data1.

\(f = StinemanInterpolation[data1]; 

    g = ShapedInterpolation[data1, Method -> Shaped]; \)

Here is the symbolic form of the shaped interpolant.

g[x]

And now we can compare the derivatives.

\({fder, gder} = D[{f[x], g[x]}, \ x]; 
    
{fder2, gder2} = D[{f[x], g[x]}, \ {x, 2}]; \)

We will get nicer plots if we do not plot through the knots. For this we introduce a function that plots on disjoint intervals.

\(DiscontinuousPlot[f_, \ {x_, \ a_, \ b_, \ c__}, \ opts___]\  := \ 
      \((
		pl[{aa_, \ bb_}]\  := \ 
          Plot[f, {x, aa + 10\^\(-6\), \ bb\  - \ 10\^\(-6\)}, 

            \ \ \ \ \ \ \ \ \ \ \ DisplayFunction\  -> \ Identity, \ 
            opts]; 
		
        pr\  = \ \(PlotRange\  /. \ {opts}\)\  /. \ 
            PlotRange\  -> \ Automatic; 
\ \ \ \ 
        Show[pl\ /@\ Partition[{a, \ b, \ c}, \ 2, \ 1], \ 
\ \ \ \ \ \ 
          Sequence\ @@\ 
            DeleteCases[
				{FilterOptions[Graphics, \ opts]}, \ 
              PlotRange -> _, \ ∞], 
		
          PlotRange -> If[VectorQ[pr], {{a, Last[{c}]}, \ pr}, pr]])
        \)\)

And we now plot the derivatives over the right part of the domain.

\(\(Show[
      GraphicsArray[
	
        \(DiscontinuousPlot[#, {x, 9, 11, 12, 14, 15}, 
			\ 
              Frame -> True, Axes -> None, \ 
			
              DisplayFunction\  -> Identity, 
		\ \ \ 
              PlotPoints -> 10, PlotRange -> {\(-5\), 60}]&\)\ /@

          \ \ \ {fder, \ gder}]]; \)\)

We can see that the first derivatives exist and are continuous in both cases, but that it undergoes very rapid change in the Stineman case (though at the knot at x = 14, the shaped derivative has a visible kink too). We can see the distinction another way by looking at the second derivatives. The values for the Stineman method get too large and also take large jumps.

\(\(Show[
      GraphicsArray[
		
        \(DiscontinuousPlot[#, {x, 6, 8, 9, 11, 12, 14, 15}, 
				
              Frame -> True, Axes -> None, 
				
              DisplayFunction -> Identity, PlotPoints -> 10, 

              			\ \ PlotRange -> {\(-130\), 160}]&\)\ /@
          \ \ {fder2, gder2}]]; \)\)

Recall that cubic splines have the property that, among all curves through a given point, they minimize the integral of the second derivative squared; but this is accomplished at the cost of the failure to preserve monotonicity. The pure second derivative is not the only measure of “fairness” of a curve (see [Roulier 1994] for a survey of measures of fairness); for example, from a physical perspective it makes sense to consider the strain energy, which is the integral of the square of the curvature, assuming the curve is reparametrized according to arc length. We will be content with just showing how the symbolic form allows us to integrating the square of the second derivative, thus getting a rough measure of the fairness of each curve. 

The package contains a special function to numerically integrate Which statements. First we must put the integrand into a Which form.

\(fder2squared = 
      Which@@\(Join@@
          \((Partition[List@@fder2, 2, 2] /. {a_, e_} 
                :> {a, e\^2})\)\); 
    gder2squared = 
      Which@@\(Join@@
          \((Partition[List@@gder2, 2, 2] /. {a_, e_} 
                :> {a, e\^2})\)\); 


    NIntegrateWhich[fder2squared, {x, 1, 15}]\)

NIntegrateWhich[gder2squared, {x, 1, 15}]

Because of figure , the comparison yields no surprise: the shaped interpolation has a smaller squared second derivative, on average, and so can be viewed as being less sharply curved.

Acknowledgements.  I am grateful to Brenda Kroschel (née Anderson) and John Roulier for telling me about many of these modern methods. Roulier has written many papers on the subject of monotonic interpolation by functions and parametric functions such as splines; see, for example, [Roulier 1994].

References

B. K. Anderson, Shape-preserving piecewise rational interpolants, masters thesis, Clemson University, 1987.
R. Delbourgo and J. A. Gregory, The determination of derivative parameters for a monotonic rational quadratic interpolant, IMA Journal of Numerical Analysis 5 (1985) 397–406.
J. Roulier, Measures of fairness for curvs and surfaces, in Designing Fair Curves and Surfaces, N. Sapidis, ed., SIAM, Philadelphia, 1994, 75–122.
R. Stineman, A consistently well-behaved method of interpolation, Creative Computing (July, 1980) 54–57.
S. Wagon, Monotonic interpolation, Mathematica in Education and Research 4:3 (1995) 49–52.

About the Author

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including the just-released second edition of Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

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