March 7, 2017 --- Class 17 --- Numerical Integration in One Dimension
Error Estimation, Simpson's Rule Rederived
Numerical Integration
We talked about some of the standard methods for numerical
integration in one dimension. This is discussed in Section 11.1
of CSM. Additional material is available in Numerical Recipes
or in Abramowitz and Stegun. We discussed the rectangle rule,
trapezoid rule, and Simpson's rule.
Simpson's rule, with Mathematica
There is a Mathematica notebook showing how integrals needed for
the derivation of Simpson's rule can be easily done. Of course,
the integrals can also be done by hand. The notebook can be found
in ~sg/Documents/SimpsonsRule2016.nb and there is a pdf in
~sg/Documents/SimpsonsRule2016.pdf.
Error estimates
We derived error estimates by using Taylor expansions of the
function, for the rectange and trapezoid rules. For the rectangle rule,
total error = (1/2) h (b-a)
where h is the bin width, sometimes called \Delta x, the integral
goes from a to b and stands for the average of the derivative
of the function.
For the trapeziod rule, the corresponding formula is
total error = (1/6) h^2 (b-a)
An interesting point about the trapeziod rule is the error only
has even terms in h, i.e., the next order term is h^4, not h^3. This
will be important later.
The derivation for Simpson's rule is a litte different in that we
just Taylor expand the function around the midpoint of the pair of
bins and get an expression that involves the second derivative of f
and indicated an error of order h^5 (in the pair of bins). Then, we
needed to use a finite difference approximation for the second
derivative. We rederived Simpson's rule and found the local error
is order h^5, including the error in the approximation to the 2nd
derivative. The global error is then order h^4. This is two orders
better than the trapezoid rule. This remarkable benefit comes from
the symmetry that arises in using two adjacent bins. Another way to
see this comes from considering what we might be neglecting when we
used a parabolic function to approximate f(x), which was the first
way we derived Simpson's rule. Our parabolic approximation to f(x)
can differ from the true function by cubic, quartic and higher order
polynomials, but all those polynomials must vanish at x0, x1 and x2.
It is convenient to translate the x axis by x1, so our integration
region goes from -h to +h. In that case, the higher order
polynomials vanish a -h, 0, and +h. The lowest order polynomial
that does so is (x+h)x(x-h). We can add another power of x to
get a higher order polynomial that vanishes at the correct places.
The Mathematica notebook mentioned above, integrates both these
functions. The integral of the cubic is zero by symmetry. The
integral of the quadratic does not vanish, but its integral is
proportional to h^5, thus confirming the local error estimate above.
Simpson's Rule Rederived
For the following discussion, let N be the number of bins.
The leading error in the trapeziod rule is 1/N^2. We can use an
idea similar to what we did to remove sample size bias, but taking
into account the fact that dependence of the error on N is different.
Starting from the trapezoid rule, we used the formulae below
for correcting the leading 1/N^2 errors.
S_N = A + E/N^2
S_2N = A + E/(2N)^2
so A = 4/3 S_2N -1/3 S_N
This technique can be extended in the method known as
Romberg integration. This expanded upon in Numerical Recipes.
Practical codes may be found there and are discussed below.
Romberg method will be detailed in the next class.
Applying the formula above to a one-dimensional integral, we found that
the weighting of the sum of two expressions on the right hand side
exactly coincides with the Simpson's rule formula. We previously
found that Simpson's Rule has a leading error of order 1/N^4. Since
the approach above removes the leading 1/N^2 error and since the
error of the trapezoid rule only has even powers of 1/N (asserted, but
not proved), this is not surprising.