<< Graphics`;<<Miscellaneous`RealOnly`

Error approximation in Taylor Series

This lab is an exploration of some of the problems that come up when you approximate a function with a finite Taylor polynomial.

WARNING:  THIS LAB IS NOT COMPLETE YET.  I AM PUTTING IT UP FOR COMMENTS AND FOR MY OWN AMUSEMENT.  I WILL EVENTUALLY POLISH THIS UP AND FINISH THE COMMENTARY, BUT THIS IS ALL THERE IS FOR NOW.  YOU HAVE BEEN WARNED...

Find the size of the interval, given the number of terms

In this section, we will examine the situation where we want to approximate our function with a given degree Taylor polynomial.  We thus need to find out what domain this polynomial is guaranteed to fall within our error tolerence.

A bounded set of derivatives (the easy case)

Let's start out with an easy case, where we know that the derivatives are all bounded (and we know that bound):

f[x_] := Cos[x]; c = π/4; nDegree = 12; 
    yRange = {-6, 6};

c is the center around which the Taylor series is expanded, nDegree is the degree Taylor function we want to use, and yRange is just the range of y values we want to see when we graph these (not used for anything else).

We want to know how large the domain of this Taylor series can be if we use a (12^th) degree Taylor approximation.

In order for the maximum error bound to be correct, you must enter a number or a function of n here:

maxNthDerivValue[n_] := 1

This is actually the most difficult part of the problem in the general case.  For cos x, of course, this is trivial (why?), but for a generic function, finding the maximum size of its (n^th) derivative is very hard (especially if you don't know in advance how many terms you will need to use).

We need to define a few functions for working with errors.  The first one just computes the actual error between the original function and the (n^th) degree Taylor Polynomial (as a function of the degree).

errorFunction[n_] := 
      N[f[x] - Normal[Series[f[x], {x, c, n}]], 20]

So, for example, in our case:

errorFunction[nDegree]

Now, we need to define a function to approximate our maximum possible error ((((M(x - c))^(n + 1))/((n + 1)!)), where M = max[(f^((n + 1)))(k)] ):

errorBoundFunction[n_, x_] := 
      N[((maxNthDerivValue[n + 1] (((x - c))^(n + 1)))/((n + 1)!)), 
        20]

A more accurate formula for the error is given by \(TraditionalForm\`∫\_c\%x\ \(\((x - t)\)\^n\/\(n!\)\) \(\(f\^\((n + 1)\)\)(t)\)\ dt\), which is sometimes useful (but generally much harder to work with), so we will define a second errorBoundFunction as:

errorBoundFunction2[n_, x_] := 
      Simplify[N[(∫_c^x)(((((x - t))^n)/(n!)) 
                  Evaluate[D[f[t], {t, n + 1}]]) \[DifferentialD]t, 
          20]]

The following function is designed to display an animation of two graphs side-by-side (see below for explanation).  Most of you will want to just execute this definition and move on; it isn't necessary for you to understand the details of this particular command.

animateTaylorWithErrors[maxDegree_] := 
      Table[Block[{$DisplayFunction = Identity}, {functions = 
              Plot[Evaluate[{f[x], Normal[Series[f[x], {x, c, n}]]}], 
                Evaluate[{x, xDomain[[1]], xDomain[[2]]}], 
                PlotRange -> Evaluate[yRange], 
                PlotLabel -> 
                  StringForm["f[x]= `` (Green), `` deg Taylor (Black)", 
                    f[x], n], PlotStyle -> {Green, Black}, 
                TextStyle -> {FontSize -> 7}]; 
            errors = 
              Plot[Evaluate[{errorFunction[n], errorBoundFunction[n, x]}], 
                Evaluate[{x, xDomain[[1]], xDomain[[2]]}], 
                PlotRange -> Evaluate[yRange], 
                PlotLabel -> "Actual error (Red), Max error bound \
(Yellow)", PlotStyle -> {Red, Yellow}, 
                TextStyle -> {FontSize -> 7}];}]; 
        Show[GraphicsArray[Partition[{functions, errors}, 2]]], {n, 1, 
          maxDegree, 1}]

Now, let's look at some graphs.  These graphs show the original function (green) and the (n^th) degree Taylor series (black) in the graph on the left and the actual error (f(x) - (Subscript[T, n])(x), in red) and the maximum error bound ((((M(x - c))^(n + 1))/((n + 1)!)), where M = max((f^((n + 1)))(k)) where k is any x value in the interval under consideration, in yellow).  Notice a few things:

1)  Observe how the Taylor polynomial gradually fits itself to more and more of the graph of f(x) as the degree of the Taylor polynomial is increased.

2)  Notice how the errors become small over larger and larger intervals as the degree of the Taylor polynomial is increased.

3)  Notice how the magnitude of the maximum error bound is indeed always at least as big as the actual error everywhere.  (Does it matter that sometimes the actual error and the maximum error bound have different signs?)  This means that if you can keep your maximum error bound small enough, you can be assured of the accuracy you want.

animateTaylorWithErrors[10]

So, to find our domain, we need to see when the error is either (10^(-9)) or -(10^(-9)) (above or below):

{Solve[errorBoundFunction[nDegree, x] == -(10^(-9))], 
      Solve[errorBoundFunction[nDegree, x] == (10^(-9))]}

We can safely ignore the "nonreal" solutions, so the domain is \(TraditionalForm\`x ∈ \((\(-0.365663224572150538938906353866719`20\), 
          1.9364595513670471492612007784`16.1851)\)\), including or excluding the endpoints, depending on whether you want to allow the error to actually equal (10^(-9)).  This gives us a radius of convergence of:

(1/2) (1.93645955136704714926120077847 -  -0.365663224572150538938906353866719)

Derivatives of any size (the hard case)

Now, let's look at something a little more challenging.  Let's look at the base function for a normal distribution:  f(x) = (e^(-(x^2))).  We would like to approximate this function using a (12^th) degree Taylor series, but we need to find out what domain this falls within our error tolerence.  (That way, we can do things like integrate it and use it in differential equations.)

f[x_] := (E^(-(x^2))); c = 0; nDegree = 12; 
    yRange = {-0.2, 1.2}; xDomain = {-6, 6};

I defined an xDomain here, even though we are trying to determine this.  This is for two reasons: 
1)  To draw the graphs
2)  Because there is no way to determine the maximum value of the (n^th) derivative without choosing an interval.  If you find the domain of our Taylor series extends outside this interval, you should make the xDomain larger and recompute the domain.

The problem here is that the derivatives are not bounded over our interval (or much of any interval, actually).  Let's examine the first 10 derivatives over our domain:

Table[
      Plot[Evaluate[D[f[x], {x, k}]], 
        Evaluate[{x, xDomain[[1]], xDomain[[2]]}], 
        PlotRange -> All, 
        PlotLabel -> StringForm["`` derivative of f(x)", k]], {k, 0, 
        20}]

Thus the maxNthDerivValue gets huge in this case (and depends on the value of n).  Ugh.

The following function attempts to use the Mathematica FindMinimum function to find the maximum value of the absolute value of the (n^th) derivative of f(x).  Unfortunately, in the most general case, this cannot be relied on to find the absolute maximum all the time (I built in a few tricks to try to find it more often).  It is important that if you use this to compute an error bound, you should check the graph of the (n^th) derivative that you decide to use to make sure you haven't missed anything important.  (It is not crucial for you to understand the details of this function.)

maxNthDerivValue[
        n_] := {Clear[x]; 
          Max[Abs[{FindMinimum[
                    Evaluate[
                      D[f[x], {x, 
                          n}]], {x, (xDomain[[2]] - xDomain[[1]])/
                          2 + xDomain[[1]], xDomain[[1]], 
                      xDomain[[2]]}, MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[D[f[x], {x, n}]], {x, xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[D[f[x], {x, n}]], {x, xDomain[[2]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[
                          f[x], {x, 
                            n}]], {x, (xDomain[[2]] - 
                              xDomain[[1]])/2 + xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[f[x], {x, n}]], {x, xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[f[x], {x, n}]], {x, xDomain[[2]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], 
                Evaluate[D[f[x], {x, n}]] /. x -> xDomain[[1]], 
                Evaluate[D[f[x], {x, n}]] /. 
                  x -> xDomain[[2]]}]]}[[1]]

maxNthDerivValue[20]

We can now graph the function, Taylor series, and errors for different values of n.  (Ignore the error messages.)

animateTaylorWithErrors[20]

Why does it seem that the Taylor series here only changes every (2^nd) degree (i.e., 3 and 4 look the same, 5 and 6 look the same, etc.)?  Are they really the same, or just "close"?

In order to find the domain that keeps the error less than (10^(-9)), we must solve:

{Solve[errorBoundFunction[nDegree, x] == -(10^(-9))], 
      Solve[errorBoundFunction[nDegree, x] == (10^(-9))]}

Recall what I said above about the necessity of checking that maxNthDerivValue is accurate for this degree Taylor series.  Let's compare the maximum value of the absolute value of the derivative we are using here (nDegree + 1) with a graph of that derivative:

maxNthDerivValue[nDegree + 1]

Plot[
      Evaluate[D[f[x], {x, nDegree + 1}]], {x, xDomain[[1]], 
        xDomain[[2]]}, PlotRange -> All]

Did it work?  What is the domain and radius of convergence?

Answer

domain = {-0.36378187103439624`, 0.36378187103439624`}

((domain[[2]] - domain[[1]])/2)

Find the number of terms, given the size of the interval

In this section, we will examine the situation where we want to approximate our function over a given interval.  We thus need to find out how many terms we need to take to guarantee our error is small enough.

A bounded set of derivatives (the easy case)

Let's redefine our function back to the cosine:

f[x_] := Cos[x]; c = π/4; xDomain = {-2.  π, 2.  π}; 
    yRange = {-6, 6};

maxNthDerivValue[n_] := 1

Notice that we specify the xDomain over which we want to use the Taylor series, but not the degree of the Taylor series (we are trying to find that).

We will make use of the errorFunction and errorBoundFunction defined above.

So, what degree Taylor series do we need to take to achieve a given level of accuracy over an interval?  If we want an accuracy of at least (10^(-9)) over our entire domain, for example, we would need to ensure that the maximum error is less than (10^(-9)) at both endpoints (why there?).  So, we can take our errorBoundFunction and substitute in the xvalue at the endpoint that maximizes (x - c).  We then have an equation in n that we can set equal to (10^(-9)) and solve.  Unfortunately, Mathematica can't solve this equation directly, so we either have to use Newton's method to solve (use FindRoot, though it may take some fooling around to find an initial guess that works) or just display a table of values and solve by inspection.

First, let's choose the x value such that (x - c) is a maximum size (why?):

Clear[n, errorBound, endPoint]; 
    endPoint = Max[Abs[xDomain - c]]

Set it equal to (10^(-9)) and solve for n:

FindRoot[errorBoundFunction[n, endPoint] == (10^(-9)), {n, 5}, 
      WorkingPrecision -> 30, MaxIterations -> 30]

WARNING:  To trust FindRoot for a calculation like this, you will probably have to increase the WorkingPrecision to be significantly more decimals than your desired error (30 digits in this case) and you may also need to set MaxIterations to a larger number than the default of 15 (if you get an error message saying "Newton's method failed to converge to the prescribed accuracy after 15 
iterations.").  This isn't a big deal, but it is an extra thing to worry about.  You could just list the errors for different values of n and find your answer by inspection:

TableForm[Table[{n, errorBound[n]}, {n, 1, 40, 1}], 
      TableHeadings -> {None, {"n", "Maximum error bound"}}]

So, taking n = (30^th) degree should do the trick.  Recall that this is the worst your error should be on the interval.  The actual error may be significantly smaller.  It is interesting to look at how the error varies with the degree of the Taylor series:

Plot[errorBound[n], {n, 0, 30}, PlotRange -> All]

Notice an interesting fact:  The error actually gets worse as you increase the degree between 1 and 5, then improves steadily.

Explain what is going on here.  How can the error be worse with more terms?  Does this mean you would be better off using a first degree series than using a (5^th) degree?  Explain why/why not.

Derivatives of any size (the hard case)

Now, let's look back at  f(x) = (e^(-(x^2))).  We would like to approximate this function using a Taylor series over a particular interval.

f[x_] := (E^(-(x^2))); c = 0; xDomain = {-5. , 5. }; 
    yRange = {-0.2, 1.2};

We also need to go back to our complicated function to find the maximum size of the (n^th) derivative.

maxNthDerivValue[
        n_] := {Clear[x]; 
          Max[Abs[{FindMinimum[
                    Evaluate[
                      D[f[x], {x, 
                          n}]], {x, (xDomain[[2]] - xDomain[[1]])/
                          2 + xDomain[[1]], xDomain[[1]], 
                      xDomain[[2]]}, MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[D[f[x], {x, n}]], {x, xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[D[f[x], {x, n}]], {x, xDomain[[2]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[
                          f[x], {x, 
                            n}]], {x, (xDomain[[2]] - 
                              xDomain[[1]])/2 + xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[f[x], {x, n}]], {x, xDomain[[1]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], FindMinimum[
                    Evaluate[-D[f[x], {x, n}]], {x, xDomain[[2]], 
                      xDomain[[1]], xDomain[[2]]}, 
                    MaxIterations -> 80, 
                    WorkingPrecision -> 30][[1]], 
                Evaluate[D[f[x], {x, n}]] /. x -> xDomain[[1]], 
                Evaluate[D[f[x], {x, n}]] /. 
                  x -> xDomain[[2]]}]]}[[1]]

Clear[n, errorBound, endPoint]; 
    endPoint = Max[Abs[xDomain - c]]

Again, let's use the errorBoundFunction, etc. defined above.

Let's take a look at the error as a function of n again (this graph takes a long time to generate on my computer, so be warned unless you have a pretty fast computer):

ListPlot[Table[{n, errorBoundFunction[n, endPoint]}, {n, 0, 250}], 
      PlotRange -> All]

(We had to use ListPlot instead of Plot because of our complicated maxNthDerivValue function.  It also caused all those error message, which are mostly harmless.)

Notice that the error bound gets extremely large before it starts to decline again for this function.  Also, all the error messages are indications where the FindMinimum function is having trouble.

In order to get our error bound within (10^(-9)), it is going to take a lot of terms.  I don't think there is a reasonable way to find this automatically, so we will just make a table of error bounds starting with the (180^th) degree and going up, finding the answer by inspection:

TableForm[
      Table[{n, errorBoundFunction[n, endPoint]}, {n, 180, 230, 1}], 
      TableHeadings -> {None, {"n", "Maximum error bound"}}]

So, it looks like it will take a (207^th) degree Taylor series to achieve this degree of accuracy throughout the domain.  Let's double-check the (208^th) derivative:

maxNthDerivValue[208]

Plot[
      Evaluate[D[f[x], {x, 208}]], {x, xDomain[[1]], xDomain[[2]]}, 
      PlotRange -> All]

Plot[Evaluate[D[f[x], {x, 208}]], {x, 4.9, 5.1}, 
      PlotRange -> All]

Did it work?

Something odd is definitely going on here.  What do you think the problem is and how would you suggest we attempt to fix it?