PROGRAMMING TIPS

Bending Plot to Your Needs

Mathematica’s new precision control allows one to generate plots that are more accurate than they deserve to be!

This column stems from two requests I received, in the same week, for advice when plotting numerically unstable functions. The examples in question were the expansion of (((1 - x))^10)and the 200th Taylor polynomial for cos x. Obtaining accurate plots of these functions is a superb exercise in the new precision control capabilities of version 3.0’s N function.

1. Difficult Polynomials

Let's start with an expanded (10^th) degree polynomial, changing the scale by (10^14).

f = Expand[(10^14) (((1 - x))^10)]

Since the unexpanded form causes no difficulty, we plot both for comparison.

\({x0, x1, y0, y1}\  = \ {0.95, 1.05, \ \(-9\), \ 9}; 

    Show[Block[{$DisplayFunction = Identity}, 

        		{Plot[f, {x, 0.95, 1.05}], 
			
          Plot[10\^14\ \((1 - x)\)\^10, {x, 0.95, 1.05}, 
				
            PlotStyle -> {AbsoluteThickness[3], \ GrayLevel[1]}]}], 
	
      PlotRange -> {\(-9\), 9}, Frame -> True, Axes -> None, 

      	Background -> GrayLevel[0.5]]; \)

The problem is subtractive cancellation caused by the very large numbers that show up as individual terms of the polynomial. Machine precision says that f(0.97) is 3.6875.

f /. x -> 0.97

But using rational arithmetic  — in essence, using infinite precision — gives us the right answer.

f /. x -> (97/100)

N[%]

We can also get the correct answer by increasing the precision of the input, 0.97, to 100 digits, for then all the arithmetic will be carried out using 100 digits.

N[f /. x -> SetPrecision[0.97, 100]]

Indeed, since we do not know in advance how many digits we might need, we can first convert the approximate real to an infinite precision number, which produces a rational number near 0.97.

approx = SetPrecision[0.97, ∞]

N[f /. x -> approx]

Note that when rationals are used, the polynomial is first evaluated using rational arithmetic throughout. This is quite slow, and we will see in a moment how to avoid it.

Now, version 3.0’s numericization function  N has some amazing precision-control capabilities (see [Sofroniou 1996]). N[x, 20] will attempt to use extra precision in the evaluation of x so that the result is guaranteed to 20 digits. This works only if the requested precision is beyond machine precision; otherwise machine precision is used throughout.

N[Sin[(10^60)], 20]

Of course, 20 digits of (10^60) would be insufficient to place the argument in a position on the number line so that its sine can be determined.

Sin[N[(10^60), 20]]

Precision[%]

Note that there is a limit on the amount of extra precision that will be used. The default setting of $MaxExtraPrecision is 50, and that would be insufficient for the following; it can be increased, but let us leave it at its default for now.

N[Sin[(10^100)], 20]

So here is how we can plot the troublesome function. We first turn the expression f into a function g, but we restrict g to operate on inexact numbers only. This means that when g sees a rational number it does nothing, thus avoining time-consuming computation with rationals. Because my Macintosh’s precision is 16, I use 17 to force the high-precision computation. $MachinePrecision will tell you the precision of your computer. And we must turn compilation off so that software arithmetic can take over.

g[x_?InexactNumberQ] := Evaluate[f]

Plot[N[g[SetPrecision[x, ∞]], 17], {x, 0.95, 1.05}, 
      Compiled -> False, Frame -> True, Axes -> None, 
      PlotRange -> {-1, 8}];

We now turn to the interesting problem of plotting high-degree, alternating Taylor polynomials. But first let us automate the process of calling on high precision. One extra twist is to make g a NumericFunction. This helps N understand that g will ultimately produce numbers and, for some complicated examples, will lead to better precision control. In fact, it is not needed for the examples presented here.

\(\(SlowButSurePlot[f_, {x_, a_, b_}, opts___] := 
      Module[{g}, 
		SetAttributes[g, NumericFunction]; 
		
        g[y_?InexactNumberQ] := Evaluate[f /. x -> y]; 
		
        Plot[N[g[SetPrecision[x, ∞]], 1 + $MachinePrecision], {x, 
            a, b}, 
			Compiled -> False, opts]]; \)\)

We first get a 200th-degree approximation to cosine. The numbers that follow show the huge variation in the summands when x is 70.

\(taylorPoly = Normal[Series[Cos[x], {x, 0, 200}]]; 

    data = Apply[List, taylorPoly] /. x -> 70.; \)

N[data
        [[{1, 2, 20, 30, 40, 50, 60, 70, -2, -1}
          ]], 3]

Using machine precision when x is 70 is not a good idea. Using 30 digits of precision is adequate for this value of x.

N[taylorPoly /. x -> {70., 70, N[70, 30]}]

Of course, a machine precision plot breaks down badly.

Plot[taylorPoly, {x, 0, 76}, PlotRange -> {-10, 10}];

But using high precision, though slow, leads to an accurate view of the polynomial.

SlowButSurePlot[taylorPoly, {x, 0, 76.1}, 
      Ticks -> {{76}, Automatic}];

The interested calculus teacher can now investigate questions such as: Where will the nth Taylor polynomial first diverge from the cosine? Of course, high precision computations are slow, and in specific cases one can often resort to some trickery to avoid them. Mathematica's symbolic power yields a closed-form expression for the Taylor polynomials of cosine.

taylorPolyFast[n_][x_]  :=  
      Evaluate[PowerExpand[
          (∑_(i = 0)^(n/2))(((((-1))^i) 
                (x^(2 i)))/((2 i)!))]]

taylorPolyFast[n][x]

The preceding expression yields an accurate plot with no numerical subtleties. Another way to avoid high precision would be to use Horner’s rule for the evaluation of the polynomial; we leave that as an exercise.

2. Overloading Plot

Another way to incorporate smart plotting is by adding a Plot option, called UseExtraPrecision, that causes a switch into high precision. While much has been written about overloading built-in functions with new cases, adding options is a little trickier. The general idea is to unprotect Plot and define it in the case that UseExtraPrecision appears as an option. In that case various modifications to the arguments should be made and then the built-in Plot gets called on the modified arguments. Of course, one modification must be to strip away the new option (the DeleteCases phrase in the code that follows), so that the second call to Plot goes through the built-in mechanism.

Unprotect[Plot]; 

UseExtraPrecision::usage = "UseExtraPrecision is an option to Plot that causes the precision of the computation to be raised if the function is numerically sensitive so that the visual results are accurate.";

Plot[f_, {x_, a_, b_}, opts___] := Block[
 {highQ, newopts, g, $MaxExtraPrecision = 200},
 highQ = UseExtraPrecision /. {opts} /. UseExtraPrecision->False;
 newopts = Sequence @@ DeleteCases[{opts}, UseExtraPrecision -> _]; 
 SetAttributes[g, NumericFunction];
 g[z_?InexactNumberQ] := Evaluate[f /. x -> z];
 If[highQ, 
   Plot[N[g[SetPrecision[x, ∞]], 1 + $MachinePrecision], {x, a, b}, 
      Compiled->False, Evaluate[newopts]], 
   Plot[f, {x, a, b}, Evaluate[newopts]]]] /; 
      MemberQ[{opts}, UseExtraPrecision->_]

Protect[Plot];

Now a call such as 

\(\(\ \ \ \ \ 
      \(Plot[Evaluate[Expand[10\^14\ \((1 - x)\)\^10]], {x, 0.95, 1.05}, 

        	\ \ \ \ \ \ \ \ \ UseExtraPrecision -> True]; \)\)\)

will yield the accurate plot shown in Figure .

This approach works, and the usage message is helpful, but we have not yet added UseExtraPrecision to Options[Plot]. Doing this is somewhat tricky because we have to make sure that the new option is not present in Options[Plot] in the cases in which it is not used, for otherwise a bad-option error will stop the plotting process. Note that the full overload requires the code that follows as well as the redefinition of Plot just given.

Unprotect[Plot];

builtInOpts = DeleteCases[Options[Plot], UseExtraPrecision -> _];

Options[Plot] = Union[builtInOpts, {UseExtraPrecision -> False}];

Plot[func_, range_, opts___] := (
  Unprotect[Plot]; oldOpts = Options[Plot];
  Options[Plot] = builtInOpts;
  result = Plot[func, range, opts];
  Options[Plot] = oldOpts; Protect[Plot];
  result) /;  !MemberQ[{opts}, UseExtraPrecision -> _]  &&
     MemberQ[Options[Plot], UseExtraPrecision -> _]
     
Protect[Plot];

Now our new option appears in the options for Plot.

Options[Plot, UseExtraPrecision]

The electronic version of the plot overload includes some additional enhancements to Plot. The ShowPlotPoints and PointStyle options show the points that Plot uses in attempting to make a smooth curve. And a case is included so that one can plot discontinuous function by using an iterator of the form {x, a, b, c, d}. An example based on the gamma function follows.

\(\(Plot[\(x!\), {x, \(-6\), \(-5\), \(-4\), \(-3\), \(-2\), \(-1\), 4}, 
      
	PlotRange -> {\(-7\), 13}, Frame -> True, 
      ShowPlotPoints -> True, 
	PlotDivision\  -> \ 3, 
	
      PointStyle -> AbsolutePointSize[3], 
	
      FrameTicks\  -> 
        \ {Range[\(-6\), 4], {\(-5\), 0, 5, 10}, \ None, None}, 
	
      GridLines 
        -> {
          \({#, {GrayLevel[0.5]}}&\)\ /@
            \ {\(-6\), \(-5\), \(-4\), \(-3\), \(-2\), \(-1\)}, \ None}, 
	
      Axes -> {Automatic, None}, \ AxesStyle -> GrayLevel[0.6]]; 
    \)\)

And we conclude with an example of a function that, I think, does not admit any tricks to avoid the numerical instability (aside from using series, but then the user would then need to work out how many terms to use). Thus some high precision approach seems necessary.

\(f = Sin[x]\ 10\^x\ Log[1 + 10\^\(-x\)]; 

    Show[Block[{$DisplayFunction = Identity}, {
			
          Plot[f, {x, 15, 20}, 
            PlotStyle -> {AbsoluteThickness[3], GrayLevel[0.6]}], 

          			Plot[f, {x, 15, 20}, UseExtraPrecision -> True]}], 
	
      Frame -> True, Axes -> None]; \)

Overloading a built-in function should probably be avoided because we really don't want Plot working differently on different machines. But now that we know how to do it, it is hard to resist! Here is one more: adding a CountEvaluations option to NIntegrate that causes the number of times the function is evaluated to be returned. The implementation is simple (though we leave as an exercise the necessary steps to add the new option to Options[NIntegrate]).

Unprotect[NIntegrate];

NIntegrate[f_, {x_, a_, b_}, opts___] := Module[{countQ, c},
  countQ = CountEvaluations /. {opts} /. CountEvaluations -> False;
  newopts = Sequence @@ DeleteCases[{opts}, CountEvaluations -> _ ];
  If[countQ, c = 0; 
     {NIntegrate[c++; f, {x, a, b}, Evaluate[newopts]], c},
      NIntegrate[f, {x, a, b}, Evaluate[newopts]]]]  /; 
                    MemberQ[{opts}, CountEvaluations -> _]

Protect[NIntegrate];

NIntegrate[Sin[Cos[x]], {x, 0, 5}, CountEvaluations -> True]

NIntegrate[Sin[Cos[x]], {x, 0, 5}, CountEvaluations -> True, 
      WorkingPrecision -> 40, PrecisionGoal -> 30]

Calculus and numerical analysis students find the evaluation counts intriguing because of the very small number of times — like 405 — that the function gets evaluated. The traditional error bound for Simpson’s rule would predict that getting an error below (10^(-30)) would require about 100,000,000 f-evaluations! The dramatic improvement is due to the adaptive algorithms underlying NIntegrate, which uses an adaptive Gaussian quadrature. In fact, even an adaptive Simpson’s rule would manage with only a few hundred f-evaluations.

Acknowledgement.  I am grateful to Michael Trott of Wolfram Research, Inc., for much valuable help with all aspects of this study.

Reference

[Sofroniou 1996]  M. Sofroniou, Numerics in Mathematica 3.0, The Mathematica Journal, 6:4, 64-73.

About the Author

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 author or
coauthor
of several books, including Mathematica in Action, Animating
Calculus, VisualDSolve, and Which Way Did the Bicycle
Go?. He also teaches a course at  Rocky Mountain Mathematica,
a high-altitude workshop held in July in the mountains of Colorado.

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

Plot Overload for Electronic Version

Needs["Utilities`FilterOptions`"];

Unprotect[Plot]; 
Clear[Plot];

Plot::usage = "Plot[f, {x, xmin, xmax}] generates a plot of f as a function of x from xmin to xmax. Plot[{f1, f2, ... }, {x, xmin, xmax}] plots several functions fi. Plot[f, {x, x0, x1, x2, ...}] superimposes plots from x0 to x1, from x1 to x2, from x2 to x3, and so on.";

ShowPlotPoints::usage = "ShowPlotPoints is an option to Plot that causes the points defining the plot to be superimposed."; 

PointStyle::usage = "PointStyle is an option to Plot and ListPlot that sets the style of the points."; 

UseExtraPrecision::usage = "UseExtraPrecision is an option to Plot that causes the precision of the computation to be raised if the function is numerically sensitive so that the visual results are accurate.";

builtInOpts = DeleteCases[Options[Plot],
     ShowPlotPoints | UseExtraPrecision | PointStyle ->_];

Options[Plot] = Union[builtInOpts,
  {PointStyle -> {PointSize[0.015]}, 
   ShowPlotPoints -> False, UseExtraPrecision -> False}];
           
Plot[f_, {x_, a_, b_}, opts___] := 
   Module[{g,sty,ptsQ, highQ,newopts, df},
    {sty, ptsQ, highQ} = 
  {PointStyle, ShowPlotPoints, UseExtraPrecision} /. 
    {opts} /. Options[Plot]; 

SetAttributes[g, NumericFunction];

g[z_?InexactNumberQ] := Evaluate[f /. x -> z];

If[!ListQ[sty], sty = {sty}]; 

newopts = Sequence @@ Join[DeleteCases[{opts}, 
  UseExtraPrecision | ShowPlotPoints | 
   PointStyle -> _], builtInOpts]; 

p = Block[{$DisplayFunction = Identity},
 If[highQ,Plot[N[g[SetPrecision[y, Infinity]],
       1 + $MachinePrecision], {y, a, b}, 
      Compiled -> False,Evaluate[newopts]], 
  Plot[f, {x, a, b}, Evaluate[newopts]]]]; 

Show[p,If[ptsQ, Graphics[Join[sty, 
  Cases[p, _Line, Infinity] /. Line[z_] :> Point /@ z], 
  FilterOptions[Graphics, newopts]],{}]]] /; 
     !FreeQ[First /@ {opts}, ShowPlotPoints|
     UseExtraPrecision | PointStyle]

Plot[f_, {x_, z__}, opts___] := Show[
  Block[{$DisplayFunction = Identity},
    Plot[f, Evaluate @ Prepend[#,x], opts]& /@ 
     Partition[{z}, 2, 1] ] ] /; Length[{z}] > 2
		
Plot[func_, range_, opts___] :=
 ( 
 Unprotect[Plot];
   oldOpts = Options[Plot];
   Options[Plot] = builtInOpts;
   res = Plot[func, range, opts];
   Options[Plot] = oldOpts;
   Protect[Plot];
   res) /;  FreeQ[ First /@ {opts}, ShowPlotPoints|UseExtraPrecision|PointStyle]  &&
     !FreeQ[Options[Plot], ShowPlotPoints|UseExtraPrecision|PointStyle]
		
Protect[Plot];