Programming Tips

Mathematica in Education and Research
Volume 4 Number 2
Copyright 1995 TELOS/Springer-Verlag

Taylor Polynomials, Efficiently

How to efficiently show a function and its Taylor polynomials.

by Stan Wagon

A common objective for computer-literate calculus teachers is to show a function together with one or more of its Taylor polynomials. While accomplishing this in particular cases is straughtforward, writing an efficient general-purpose routine brings in several subtle points about Mathematica and symbolic computation. The programming ideas discussed here arose in connection with the Taylor polynomial routine in one of the labs in [PW].

Of course, the first order of business is to get the Taylor series; not surprisingly Series is the thing to use.

Series[Sin[x], {x, 0, 7}]

Mathematica keeps track of series in a special form (FullForm[Series[]]; use ?SeriesData to get more information about this internal issue); to get rid of the error term and get a polynomial, just wrap Normal around the output.

Normal[%]

Now, we will want all the partial sums of this polynomial and we wish to avoid using Series again, since all the information we need is already computed. Since the polynomial is a Plus object (its internal form is Plus[ , , , . . .]), we can apply List to it; this will replace the Plus with List. FoldList then gets us all the partial sums. Note that this naturally avoids including any terms with coefficient 0, such as the even powers of x in this case.

Apply[List, %]

FoldList[Plus, First[%], Rest[%]]

Plot[Evaluate[Prepend[%, Sin[x]]], {x, - Pi, Pi}];

Figure 1.  The sine function together with some of its Taylor polynomials.

Now, there are at least three complications that can arise as we try to turn this basic idea to a general routine.

-The routine should work even when the input is simply a monomial such as 17x2. Our first approach fails because the series would not be a Plus-object.

- The routine should include the degree-0 approximation even if that happens to be zero. As given above, such a zero term would not appear.

- We should allow the possibility of an animation that shows one polynomial at a time, with a label telling the degree of the polynomial.

The monomial issue is not too hard to deal with. If ser denotes the series, then we must distinguish between ser being of the form Plus[], or of another form (Times[] or Power[]). In the latter case, we want just {ser}. We can therefore proceed by checking whether Head[ser] is or is not "identical to" (===) Plus. And we can take care of the zero-constant-term issue by prepending f[0] and using Union to eliminate the duplication that this can cause.

Clear[f, TaylorPolys];
TaylorPolys[f_, {x_, c_, nmax_}] := (
  ser =  Normal[Series[f, {x, c, nmax}]];
  P = If[Head[ser] === Plus, Apply[List, ser], {ser}];
  P = Union[Prepend[P, f /. x->c]];
  FoldList[Plus, First[P], Rest[P]])

TaylorPolys[Sin[x], {x, 0, 5}]

TaylorPolys[17 x^2, {x, 0, 5}]

But there are still serious problems. If f is Sin[x] + Cos[x] and the expansion is centered at 1, then the routine fails miserably because the constant term is Sin[1] + Cos[1] and, in the sorting of the output, this gets split up with Sin[1] being placed at the end; our code above assumed that the constant term was the lead-off term.

TaylorPolys[Sin[x] + Cos[x], {x, 1, 1}]

Normal[Series[Sin[x]+Cos[x], {x, 1, 1}]]

Fortunately, there is an easy fix (this is often the case ± finding the glitches in the first place is the hard part!). As with so many programs of this type, the programmer must protect the code form users (including the programmer!) who neglect to use numerical arguments. Simply using N[] at a couple of key points provides such protection.

TaylorPolys[f_, {x_, c_, nmax_}] := Module[{cc = N[c], ser, P},
  ser =  Normal[Series[N[f], {x, cc, nmax}]];
  P = If[Head[ser] === Plus, Apply[List, ser], {ser}];
  P = Union[Prepend[P, N[f] /. x->cc]];
  FoldList[Plus, First[P], Rest[P]]]
  
TaylorPolys[Sin[x] + Cos[x], {x, 1, 1}]

The final subtlety concerns the determining of the degree of a polynomial. Exponent does the job, though beware: Exponent[2 + (x - 2)^3, x - 2] does not work, but Exponent[2 + (x - 2)^3, x] does return the proper answer, 3. And since Exponent[0, x] returns -Infinity when we would prefer to consider the degree of 0 as being 0, we follow up with the substitution -Infinity -> 0.

The complete and, I hope, bug-free, code for generating either a still showing a function with one of its Taylor polynomials, or a movie showing several of the Taylor polynomials one at a time is given at the end of the column. A natural extension of this code would be to generate a single image with several Taylor polynomials; we leave this as an exercise. Note the use of StringForm to make a label that incorporates the degree of the polynomial. And note how Scan is used: we first define and name the plot of f; then we scan the list of Taylor polynomials, showing each one together with the plot of f. And finally, it is important when generating animations to control the plot ranges: they should be the same in all frames. We accomplish that by first getting the range for the plot of the function alone (via 
PlotRange[Plot[f, DisplayFunction->Identity,...]]; the Identity setting suppresses the display here). We then use these numbers when we actually show the images later.

Here are some sample examples. The first one shows a familiar case; convergence happens over all of (—®, ®) of course.

TaylorMovie[Sin[x], {x, -6, 6, 0}, 5];

Figure 2.  The sine function and its 5th-degree Taylor polynomial.

The case that follows shows that our routine is adept at dealing with symbolic input.

TaylorMovie[Sin[x] + (2/Pi) Cos[2 x],
	{x, -10, 10, 1}, 11];

Figure 3.  The degree-11 approximation to Sin[x] + (2/Pi) Cos[2 x].

Figure 4 again shows a typical case for which the convergence is over all of  (—®, ®).

TaylorMovie[Exp[-x^2], {x, -2, 2, 0}, 4]
TaylorMovie[Exp[-x^2], {x, -2, 2, 0}, 16]

Figure 4.  The degree-4 and degree-16 approximations to Exp[-x^2] lead one to think, correctly in this case, that the convergence occurs over the whole real line.

Examining several functions such as exp(-x2) might lead students to think that most continuous functions have Taylor polynomials that converge throughout their domain. Of course, this is false; the arctangent function is a convenient example, since there is no convergence outside of (-1, 1). Examining several Taylor approximations (see Figure 4) shows that convergence seems to be prevented from passing through an invisible barrier at |x| =1. I say "invisible" because the graph of the arctangent function gives no clue that such a barrier should exist. Complex numbers explain the situation, because the derivative of arctan(x) is 1/(1 + x2), which, in the complex plane, has a singularity at –i; this prevents the complex circle of convergence from having radius 1.

TaylorMovie[ArcTan[x], {x, -3, 3, 0}, 3]
TaylorMovie[ArcTan[x], {x, -3, 3, 0}, 15]

Figure 4.  Two Taylor approximations to the arctangent function show that a barrier seems to block convergence outside of the interval (—1, 1). This barrier is really the ghost of the imaginary number i making its presence felt, because the derivative of the arctangent function is undefined at –i.

Exercise.   Modify the TaylorMovie code so that it allows the user the option of superimposing several Taylor polynomials on one graph.

TaylorMovie code

TaylorMovie::usage = "TaylorMovie[f, {x, a, b, c}, {nmin, nmax}, opts] \
generates images of the Taylor polynomials for f at x = c between a and \
b. The first and last degrees are given by nmin and nmax. Regular \
Graphics options may be added. For a single image use simply: 
TaylorMovie[f, {x, a, b, c}, n]."; 

TaylorMovie[f_, {x_, a_, b_, c_}, {nmin_, nmax_}, opts___] := Module[
 {P, plRange, funcPlot, y, cc = N[c], ser},

 ser = Normal[Series[N[f], {x, cc, nmax}]];
 P = If[Head[ser] === Plus, List @@ ser, {ser}];
 P = Union[Prepend[P, N[f] /. x->cc]];
 P = FoldList[Plus, First[P], Rest[P]];

 plRange=PlotRange[funcPlot = Plot[f, {x, a, b},
   DisplayFunction->Identity, PlotStyle->AbsoluteThickness[0.5]]];

 Scan[(
   polyDegree = Exponent[#, x] /. -Infinity -> 0;

   If[nmin <= polyDegree <= nmax || 
   (nmin == nmax && polyDegree[#] == polyDegree[Last[P]]), 

   Show[funcPlot,
	 Plot[#, {x, a, b},
		DisplayFunction->Identity, PlotPoints->100,
		PlotStyle->AbsoluteThickness[1.5]],
	 opts,
	 Epilog->{AbsolutePointSize[7], Point[{cc, f /. x->cc}]},
	 PlotRange->plRange,
 	PlotLabel->
	   FontForm[
	      StringForm["Taylor polynomial of degree ``", polyDegree],
	   {"Times-Bold",10}],
	 DisplayFunction->$DisplayFunction]]) &,  P]]
 
 TaylorMovie[f_, {x_, a_, b_, c_}, n_, opts___] :=
     TaylorMovie[f, {x, a, b, c}, {n, n}, opts]

References

Ed Packel and Stan Wagon, Animating Calculus: Mathematica Notebooks for the Laboratory, W.˚H. Freeman, New York, 1994.