Exploring errors in approximating definite integrals numerically

This lab is designed to help you investigate how accurate the different standard methods of numerically approximating a definite integral are and how "efficient" each of these is (informally).  

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...

The function and its integral (run this first)

For the purposes of this lab, I have chosen an integral that can be integrated exactly.  This allows us to actually compare the computed values to the true value of the definite integral.  (I did try to pick one that most of us would rather not try to actually integrate, however...)  

First, some preliminary commands to set things up.  The first one tells Mathematica to only use Real numbers (i.e., not complex numbers).

<< Miscellaneous`RealOnly`

This loads the basic graphics libraries (I generally load this one "just in case" I need it):

<< Graphics`

Now, we define our function.  If you want to try everything again with a different function, this is where you need to make that change.

f[x_] := (Sqrt[Sin[(x^2)] + 1])

We want to integrate (∫_a^b) f(x) dx, so let's define our limits of integration:

a = 0; b = π/2;

Let's graph it.  (Warning:  If you do change the function,you will probably want to change the PlotRange as well.)

Plot[f[x] , {x, a, b}, PlotRange -> {0, 1.5}]

Okay, now we want to see what the indefinite integral is (assuming it can actually be found).

indefIntegral = ∫f[x] \[DifferentialD]x // 
        FullSimplify

So, the actual value of the definite integral can be found to 20 digits accuracy with the following command (note that we give this a name so we can refer back to it later).  Notice that I use the N[...,20] command repeatedly below to specify that I want Mathematica to compute everything with 20 digits of precision.

actualValue = 
      N[(∫_a^b) f[x] \[DifferentialD]x, 20]

The left and right-hand sums

Defining the Left/Right sums

The simplest way to approximate an integral is using the left-hand sum or the right-hand sum.

leftSum[n_] := {Δx = ((b - a)/n); 
        N[(∑_(k = 0)^(n - 1))(f[
                  a + k Δx] Δx), 
          20]}

rightSum[n_] := {Δx = ((b - a)/n); 
        N[(∑_(k = 1)^n)(f[
                  a + k Δx] Δx), 
          20]}

In these definitions, I compute Δx inside the leftSum/rightSum functions because it depends on n (which you won't know until the function is called).  Also, don't confuse N with n.  N is a Mathematica function that is telling the program how many digits precision we should use (I arbitrarily chose 20 digits).  n is the number of subdivisions.

This lists the left-hand and right-hand sums for different numbers of subdivisions n:

TableForm[Table[{n, leftSum[n], rightSum[n]}, {n, 100, 1000, 100}], 
      TableHeadings -> {None, {"n", "LHS", "RHS"}}]

Error approximation for the Left/Right sums

In order to find an error bound for these approximations, it can be shown that the maximum error is given by \(TraditionalForm\`\(\(|\)\(E\_n\)\(|\)\(\(≤\)\(\ \)\(\(\(\((b - a)\)\^2\) 
                maxFirstDeriv\)\/\(2  n\)\)\)\)\), where maxFirstDeriv is the maximum size of the first derivative on the interval [a, b]. To compute this, let's make some definitions:

firstDeriv[x_] = D[f[x], x]

You can use Mathematica's FindMinimum function to find the minimum and maximum values of the derivative over the interval, but this may miss a critical point (or an endpoint) in some cases, so it is a really good idea to look at a graph of the derivative over [a, b] to check that you got the correct value.  (If you aren't too particular, you could use the trace tool to find the approximate coordinates strictly from the graph.  Click once on the graph, then move the mouse over the graph while holding down the Control key.  The coordinates of the mouse cursor are shown at the lower left edge of this window.  Not very accurate, but fast...)

Plot[firstDeriv[x], {x, a, b}, PlotRange -> All]

To get more accuracy, we need to compare 4 different values:  FindMinimum of firstDeriv, FindMinimum of -firstDeriv (Mathematica doesn't have a "FindMaximum" command, oddly enough), firstDeriv(a), and firstDeriv(b).  The largest one of these will be the maximum.  (If you know that your function f(x) is simple enough, there are snazzier ways to do this, but this works.)  The FindMinimum command I use below starts looking at x = a and looks throughout the interval [a, b].

localMin = FindMinimum[firstDeriv[x], {x, a, a, b}]

This doesn't find anything because there are no local minima to find (it doesn't look for endpoints usually).  Now let's check for a maximum:

localMax = FindMinimum[-firstDeriv[x], {x, a, a, b}]

We do get an answer here (notice that it has the wrong sign; it actually occurs at 0.512394).  Now, let's compare this with the endpoints and pick the maximum size (hence the absolute values).  The notation localMax[[1]] is used to pull just the y value out (and ignores the x value).

Max[
        Abs[{localMin[[1]], localMax[[1]], firstDeriv[a], 
            firstDeriv[b]}]] // N

Since strange things like the symbolic answer above sometimes come up, I will just copy and paste the answer to the following definition (note that you need to change this if you change the function or the interval).

maxFirstDeriv = 0.9628533646753205

This defines a function that computes the error bound for n rectangles:

maxLRError[
        n_] := (((((b - a))^2) maxFirstDeriv)/(2  n))

Compare how the error bounds ("Max error") change as you increase the number of subdivisions below:

\(TableForm[
      Table[{n, leftSum[n], maxLRError[n], 
          Abs[leftSum[n] - actualValue]}, {n, 100, 1000, 100}], 
      TableHeadings -> {None, {"n", "LHS", "Max error", "Actual error"}}]\)

The table above shows the number of subdivisions increased by an order of magnitude.  Explain how the error bound changes as a result of this.  Do you expect this to hold true for other functions and/or intervals of integration as well?  Why or why not?

How does the error bound compare to the actual error in this table?  Do you expect that relationship to be the same if you are integrating some other function?  Why or why not?  Does it decrease in the same way that the error bound does as you increase the number of subdivisions?  Does this depend on the specific function used?

We can also graph the error bound as a function of the number of subdivisions (holding everything else constant):

Plot[maxLRError[n], {n, 100, 1000}, PlotStyle -> {Red}]

What general conclusions can you draw about the accuracy and "efficiency" of the left-hand and right-hand sums in computing an integral?  For the given integral, where would you say your point of "diminishing returns" would be reached (i.e., if you increase the number of subdivisions above this point, you have to work really hard to get just a little more accuracy)?  If you stopped there, how accurate would your integral be?

The Trapezoid and Midpoint rules

Defining the Trapezoid/Midpoint rules

We can define the Trapezoid rule as follows:

trapRule[n_] := {Δx = ((b - a)/n); 
        N[Δ(x/2) (f[
                  a] + (∑_(k = 1)^(n - 1))(2 
                       f[a + k Δx]) + f[b]), 
          20]}

Define a similar function for the midpoint rule, calling it midPointRule[n].

Use a table to show the values of the Trapezoid and Midpoint rules for at least one order of magnitude of n (i.e., 100 to 1000, for example).

Error approximation for the Trapezoid/Midpoint rules

The maximum error bounds for the Trapezoid and Midpoint rules are given by the following formulae:

maxTrapError (n)=(((((b - a))^3) M)/(12  (n^2))) and \(maxMidPtError \((n)\) = \(\(\((b - a)\)\^3\) M\)\/\(24  n\^2\)\), where \(TraditionalForm\`M = 
        Max\_\(x ∈ \([a, b]\)\) | 
          f\*"""\((x)\) | \)

(Notice that these can not just be cut and pasted to define them; you must adapt them for Mathematica.)

Define two error functions (as above), maxTrapError[n] and maxMidPtError[n].  Make a table comparing these errors for different values of n (at least an order of magnitude) with the actual errors for the appropriate rules.

IF you have defined your error functions like I told you to above, you should be able to execute the following command (you may need to change the domain for n, if you used different values) to plot the two errors as a function of n:

Plot[{maxTrapError[n], maxMidPtError[n]}, {n, 100, 1000}, 
      PlotStyle -> {Green, Blue}]

How does this graph compare with the similar graph for the Left/Right sum?

Simpson's Rule

Defining Simpson's rule

Define a similar function for Simpson's rule, calling it simpsonsRule[n].

Hint:  The tricky part here is deciding how to get the coefficients to alternate 1, 4, 2, 4, ... , 2, 4, 1.  You can take care of the endpoints like we did above with the Trapezoid rule.  Also notice that 3 + (((-1))^even) = 4 and 3 + (((-1))^odd) = 2.

Use a table to show the values of Simpson's rule for at least one order of magnitude of n (i.e., 100 to 1000, for example).

Error approximation for Simpson's rule

The maximum error bound for Simpson's rule is given by the following formula:

\(maxSimpsonError \((n)\) = \(\(\((b - a)\)\^5\) M\)\/\(180  n\^4\)\), where M = 
        (Subscript[Max, x ∈ [a, b]]) | (f^((4)))(
            x) | 

(Notice that this can not just be cut and pasted to define it; you must adapt it for Mathematica.)

Define an error function (as above), maxSimpsonError[n].  Make a table comparing this error for different values of n (at least an order of magnitude) with the actual errors for Simpson's rule.

IF you have defined your error functions like I told you to above, you should be able to execute the following command (you may need to change the domain for n, if you used different values) to plot all three errors as a function of n:

Plot[{maxTrapError[n], maxMidPtError[n], maxSimpsonError[n]}, {n, 100, 
        1000}, PlotStyle -> {Green, Blue, Gold}]

Compare the relative "efficiency" of the three different approximation methods.

Computation times (the cool part)

WARNING:  If you want to actually run the examples below (as opposed to just looking at the numbers I give), you must make the definitions I ask for in the exercises above (maxMidPtError, maxSimpsonError, maxTrapError, etc.).  If things don't work below, it is probably becuase you made a mistake on one or more of those or you forgot to execute them.  You have been warned...

In order to get a better feel for how much effort it takes to use each of these methods, let's make some estimates about how long it actually takes to compute our definite integral to a desired degree of accuracy.  If we want to make our approximation accurate to at least (10^(-12)), for example, we can define:

maxErrorGoal = 1.  ((10^(-12)))

(Warning:  We defined our definitions above to have a precision of 20 digits.  If you want more accuracy than that, you will need to change those numbers in the definitions above.)

Using the left-hand (or right-hand) sum, we would need the following number of rectangles:

subdivLR = Solve[maxLRError[n] == maxErrorGoal]

Using Simpson's rule, on the other hand, we only need:

subdivSimpson = 
      Solve[maxSimpsonError[n] == maxErrorGoal]

(Why do we get negative and complex solutions here?)

This is a huge difference.  However, is this difference going to be noticeable to you if you have the computer do all the arithmetic for you?  Let's do some very unscientific experimentation to find out.  Mathematica has a function called Timing that allows you to time how much CPU time your computer needs to perform a calculation.  What we are going to do is to sample this time for different numbers of subdivisions (10 through 10010 subdivisions, increasing by 1000 each time).  Then, we will fit this data to a line so we can predict how long it would take this computer (the one you ran the test on) to compute one of these approximations for any arbitrary number of subdivisions.  (Disclaimer:  There are a lot of things that can affect this timing data, depending on what else the computer is doing while it is running.  I make no claims about any real accuracy for this.  However, it should be good enough to give us some idea of the relative effort involved in using these different approximations.  That's good enough for me...)

Let's assume all these rules take about the same amount of time to compute the same number of subdivisions.  (If you think about it, you will see that Simpson's rule isn't all that much more complicated than the left-hand rule to compute.  There are slight differences, but it takes too long on my aging home computer to get good timing data for more than one test, so let's just get some timing data for the left-hand rule.)

runTimeData = 
      Table[{n, Timing[leftSum[n]][[1]]/Second}, {n, 10, 10010, 
          1000}]

Now, we can fit this data to a line to develop a function that tells us approximately how long (in seconds) it would take to compute one of these rules for n subdivisions:

timingFunction[n_] := 
      Fit[runTimeData, {1, x}, x] /. x -> n

If you are interested, it is:

timingFunction[n]

Is it reasonable to fit this data to a line?  This makes sense theoretically (think about it), and when we plot the timing data we obtained above along with our timingFunction, it at least looks reasonable.

dataPlot = ListPlot[runTimeData, PlotStyle -> PointSize[0.02]]; 
    linearFit = Plot[timingFunction[x], {x, 0, 10000}]; 
    Show[linearFit, dataPlot]

So, to compute our integral to the required degree of accuracy using the left-sum (or right-sum), it would take us:

timingFunction[n] Second /. subdivLR

Or, if that doesn't mean much to you, it would take:

timingFunction[n]/60/60/24/365 Years /. 
      subdivLR

Simpson's rule, on the other hand, takes (ignore the negative and Nonreal answers):

timingFunction[n] Second /. subdivSimpson

Pretty impressive, no?  Now, let's take a look at time needed for computation as a function of the number of decimal places of accuracy you need (i.e., as a function of k, where (10^(-k)) is your required error bound).

This defines functions to compute the time it takes to compute integrals to (10^(-k)) accuracy for each of our methods (these formulas are kind of complicated looking because I am having Mathematica only choose the positive answers; it isn't necessary for you to understand the details of these definitions in order to understand the results):

timeLR[k_] := 
      timingFunction[
          n] /.  {n ->  Select[
                n /. Solve[maxLRError[n] == (10^(-k))], # > 
                    0 &][[1]]}

timeTrap[k_] :=  
      timingFunction[
          n] /.  {n ->  Select[
                n /. Solve[maxTrapError[n] == (10^(-k))], # > 
                    0 &][[1]]}

timeMidPt[k_] :=  
      timingFunction[
          n] /.  {n ->  Select[
                n /. Solve[maxMidPtError[n] == (10^(-k))], # > 
                    0 &][[1]]}

timeSimpson[k_] :=  
      timingFunction[
          n] /.  {n ->  Select[
                n /. Solve[maxSimpsonError[n] == (10^(-k))], # > 
                    0 &][[1]]}

Now let's look at a table of this information (you will probably need to scroll to the right to see the whole table).  Note that all the time values are in seconds here.

\(TableForm[
      Table[{10. \^\(-k\), timeLR[k], timeTrap[k], timeMidPt[k], 
          timeSimpson[k]}, {k, 1, 20}], 
      TableHeadings -> {None, {"Max error", "Left/right sum", "Trapezoid rule", "Midpoint rule", "Simpson's rule"}}]\)

How worried should you be about the negative times in the table above?  What is causing those obviously anomalous results?  (Or, I suppose it is possible that I have just invented time travel...  I doubt it, but if I did, I want full credit and a cut of the profits...)

What conclusions can you draw about the relative effieiency of these different methods?  Please try to say something a little more detailed than "Simpson's rule totally blows the others away"...

We could also plot the time necessary to compute the integral to within an accuracy of (10^(-k)) as a function of k (i.e., the number of decimal places, basically).  In order to make all the graphs visible, I am only plotting times less than one hour (red is left-hand sum, green is trapezoid rule, blue is midpoint rule, and gold is Simpson's rule).  Let's face it, if it takes more than an hour to compute, most of us will just give up and shut it down anyway.

Plot[{timeLR[k], timeTrap[k], timeMidPt[k], timeSimpson[k]}, {k, 0.1, 
        20}, PlotRange -> {0, 60. *60}, 
      PlotStyle -> {Red, Green, Blue, Gold}]

What conclusions can you draw about the relative "efficiency" of the different methods?

What is the practical limit on the accuracy achievable via each method?  (I arbitrarily decree one hour to be the longest time I am willing to wait on an answer.  You could, of course, chose a longer or shorter time.)

How accurate could you make this integral using Simpson's rule if you allow it to run for an hour?  Do you really believe this number?  Why or why not?