In numerical analysis, Quadrature is the estimation of an integral. We may want to integrate some function f(x) or a set of tabulated data. The domain might be a finite or infinite interval, it might be a rectangle or irregular shape, it might be a three dimensional cell.
We consider some simple tests of the precision of a rule, and then describe (simple versions of) the midpoint and trapezoidal rules.
We consider the Newton-Cotes and Gauss-Legendre families of rules, and discuss how to get more accurate approximations by increasing the order of the rule or subdividing the interval. Finally, we take a look at error estimation.
The "MATLAB Sum" method (I made this name up) for estimating an integral is to evaluate the function at N evenly spaced points in [A,B] and multiply by h, the spacing between the points.
Write a routine called matsum.m, which looks like
function quad = matsum ( f, a, b, n )The MATLAB Sum method might look like something like this:
h = ( b - a ) / ( n - 1 )
xvec = ?
fvec = ?
quad = h * sum ( fvec )
Exercise - Use your MATLAB sum routine to estimate the integral of 5*x4 over [0,1], using the given values of N, and record the error.
N h "MATLAB Sum" Result Error
2 1 ---------------------- ----------------------
11 0.1 ---------------------- ----------------------
101 0.01 ---------------------- ----------------------
1001 0.001 ---------------------- ----------------------
How fast is the error decreasing with h?
The MATLAB Sum rule is easy to use. It's not very accurate - it doesn't even integrate f(x)=1 exactly! We will keep this method in mind as a basis for comparison.
The midpoint rule on a single interval approximates an integral as the product of the interval width (b-a) times the function value at the midpoint:
h = ( b - a )
x = ( a + b ) / 2
quad = h * f(x)
It's a little tricky to do a composite midpoint rule with MATLAB, because we need to compute the locations of the midpoints. The easiest way to compute the midpoints is to compute the N endpoints and then average pairs of them, giving a new vector of N-1 midpoints.
Write a routine called midsum.m, which is declared as:
function quad = midsum ( f, a, b, n )The coding might look like:
h = ( b - a ) / ( n - 1 )
xvec = ?
xvec = 0.5 * ( xvec(1:n-1) + xvec(2:n) ) replaces ends by midpoints
fvec = ?
quad = h * sum ( fvec )
Exercise - Use your midpoint method to estimate the integral of 5*x4 over [0,1], using the given values of N, and record the error.
N h Midpoint Result Error
2 1 ---------------------- ----------------------
11 0.1 ---------------------- ----------------------
101 0.01 ---------------------- ----------------------
1001 0.001 ---------------------- ----------------------
How fast is the error decreasing with h?
If a quadrature rule can compute exactly the integral of any polynomial up to some specific degree, we will call this its degree of precision. Thus a rule that can correctly integrate any cubic, but not quartics, has precision 3.
To determine the degree of precision of a rule, we might look at the approximations of the integrals from 0 to 1 of the functions
p0(x) = 1since the exact value of every one of these integrals is 1.
p1(x) = 2 * x
p2(x) = 3 * x2
...
pk(x) = ( k + 1 ) * xk
Exercise - To study the precision of the midpoint method, use a single interval ( N = 2 ), and estimate the integrals of the test functions over [0,1]. The correct answer should be 1 each time.
f Midpoint Result Error
1 ---------------------- ----------------------
2 * x ---------------------- ----------------------
3 * x^2 ---------------------- ----------------------
4 * x^3 ---------------------- ----------------------
When do the results stop being essentially exact? So what is the
precision of the midpoint rule?
The trapezoidal rule over one interval multiplies the interval width times the average of the function values at the endpoints:
h = ( b - a )
Q = h * ( f(a) + f(b) )/2
To apply the composite trapezoid rule, we need to generate N points and evaluate the function there. Then we almost repeat the "MATLAB Sum" formula, except that the first and last points have a factor of (1/2) instead of 1. This amounts to a small correction to the "MATLAB Sum" formula.
Write a routine called trapsum.m, which is declared as:
function quad = trapsum ( f, a, b, n )The coding might look like something like this:
h = ( b - a ) / ( n - 1 )
xvec = ?
fvec = ?
quad = h * ( sum ( fvec ) - 0.5 * ( fvec(1) + fvec(n) ) )
Exercise - Use the trapezoid method to estimate the integral of 5*x4 over [0,1], using the given values of N, and record the error.
N h Trapezoid Result Error
2 1 ---------------------- ----------------------
11 0.1 ---------------------- ----------------------
101 0.01 ---------------------- ----------------------
1001 0.001 ---------------------- ----------------------
How fast is the error decreasing with h? Compare these results
to your MATLAB Sum results, which only differ in the "small correction".
A Newton-Cotes formula uses the idea that if it's reasonable to approximate a function by its interpolant, then it's reasonable to approximate the integral of that function with the integral of its interpolant. For instance, the trapezoidal rule is the Newton-Cotes formula of order 2, derived by integrating the linear interpolant to a function.
If our Newton-Cotes formula (over a single interval) samples the function at N points, then the interpolant will be exact for any polynomial of degree N-1 or less, hence the integral will be approximated exactly. Therefore, a Newton-Cotes rule of order N has precision at least N-1.
We didn't actually prove it can't be better. And in fact, the rules of order N have precision N when N is odd. Hence, you'd expect the midpoint rule to have precision....???
Write a routine called ncsum.m, which is declared as:
function quad = ncsum ( f, a, b, n )The coding might look like something like this:
h = ( b - a ) Watch out! This is not divided by (n-1)!
x = linspace ( a, b, n )
wvec = nc_weight ( n )
fvec = ?
quad = h * wvec' * fvec' Transpose both.
The routine nc_weight.m returns the Newton-Cotes coefficients
for the rule of given order. It is available from the web page.
Exercise - Use Newton-Cotes methods to estimate the integral of 5*x4 over [0,1], using the given values of N, and record the error.
N h Newton-Cotes Result Error
2 1.0 ---------------------- ----------------------
3 0.5 ---------------------- ----------------------
4 0.33 ---------------------- ----------------------
5 0.25 ---------------------- ----------------------
Unfortunately, "more precise" does not guarantee "more accurate"! If one rule is "more precise" than another, then it can integrate exactly polynomials of higher degree. So you might think that, for any function f(x), if we want to approximate the integral really well, we could just use Newton-Cotes methods of higher and higher precision.
But the Newton-Cotes rules do not have the very desirable property that the sequence of approximations to the integral is guaranteed to converge to the integral. It is time to examine one of our favorite examples.
Exercise: consider the Runge function 1/(1+x^2) on the interval [-4,4]. The exact integral is 2*atan(4). Use rules of increasing order, each time recording the quadrature error.
N h Newton-Cotes Result Error
3 4.0 ---------------------- ----------------------
5 2.0 ---------------------- ----------------------
9 1.0 ---------------------- ----------------------
17 0.5 ---------------------- ----------------------
Using more precise Newton-Cotes rules means using increasing degree polynomial interpolation with equally spaced nodes. When did we last try this idea, and what was the result then?
Back in the old days, arithmetic was done with a device called a "pencil", and people did everything they could to find the most effective way of reducing the pain of using this device. Gauss, in particular, had to compute many approximate integrals. He knew about the Newton-Cotes formulas, but discovered that he could approximate an integral using far fewer points (saving many, many pencils), if he was willing to replace the equally spaced sampling points by points whose location was determined in another way.
(Gauss's points are actually the places where the Legendre polynomial is zero. If you get a chance, plot the Legendre polynomial of degree 5 and then look at the values of the 5 Gauss-Legendre abscissas in the routine gl_space.m which I am about to tell you about.)
Write a routine called glsum.m, which is declared as:
function quad = glsum ( f, a, b, n )The coding might look like something like this:
h = ( b - a ) / 2 Notice this h is different!
wvec = gl_weight ( n )
xvec = gl_space ( a, b, n )
fvec = ?
quad = h * wvec' * fvec' Transpose both
The weights are computed by the routine gl_weight.m, and the
location of the points by gl_space.m, available from
the web page.
Exercise - Use Gauss-Legendre methods to estimate the integral of 5*x4 over [0,1], using the given values of N, and record the error.
N Gauss-Legendre Result Error
1 ---------------------- ----------------------
2 ---------------------- ----------------------
3 ---------------------- ----------------------
4 ---------------------- ----------------------
How fast is the error decreasing with N?
Unlike Newton-Cotes, Gauss-Legendre rules do have the important property that, in general, using the sequence of increasing-order approximations to the integral will converge to the integral.
Exercise: Vanquish the beast! Consider the Runge function 1/(1+x^2) on the interval [-4,4]. The exact integral is 2*atan(4). Use rules of increasing order, each time recording the quadrature error.
N Gauss-Legendre Result Error
1 ---------------------- ----------------------
2 ---------------------- ----------------------
4 ---------------------- ----------------------
8 ---------------------- ----------------------
As long as we keep the order low, there's nothing wrong with a Newton-Cotes rule. But if we want a more accurate approximation to an integral, we can't increase the order of the rule, because the results go bad. The solution is to use a composite rule. We break the interval [A,B] up into M subintervals and use the rule repeatedly over each subinterval. If the rule uses N points, then we want a total of M*(N-1)+1 evenly spaced points. The easiest way to organize this computation is simply to compute the endpoints of the subintervals, and call the routine we already wrote.
Write a routine, comp_ncsum.m, which looks like
function quad = comp_ncsum ( f, a, b, n, m )which approximates the integral of f(x) over [A,B] using a Newton-Cotes rule of order N, breaking the interval into M subintervals. The routine might look like this:
ends = linspace ( a, b, m+1)
quad = 0
for i = 1 to m
quad = quad + ncsum ( f, ends(i), ends(i+1), n )
Note that we are computing M+1 subinterval endpoints, which
define M subintervals.
Exercise: consider the Runge function 1/(1+x^2) on the interval [-4,4]. The exact integral is 2*atan(4). Try the composite rules of order N, with M intervals, each time recording the error.
N M Composite Newton-Cotes Error
2 4 ---------------------- ----------------------
2 8 ---------------------- ----------------------
2 16 ---------------------- ----------------------
3 4 ---------------------- ----------------------
3 8 ---------------------- ----------------------
3 16 ---------------------- ----------------------
A composite Gauss rule can be put together in much the same way as a Newton-Cotes method, except that the spacing of the points is not even.
Write a routine, comp_glsum.m, which looks like
function quad = comp_glsum ( f, a, b, n, m )which approximates the integral of f(x) over [A,B] using a Gauss-Legendre rule of order N and which breaks the interval into M subintervals. The routine might look like this:
ends = linspace ( a, b, m+1)
quad = 0
for i = 1 to m
quad = quad + glsum ( f, ends(i), ends(i+1), n )
Exercise: consider the Runge function 1/(1+x^2) on the interval [-4,4]. The exact integral is 2*atan(4). Try the composite Gauss-Legendre rules of order N, with M intervals, each time recording the quadrature error.
N M Composite Gauss-Legendre Error
2 4 ---------------------- ----------------------
2 8 ---------------------- ----------------------
2 16 ---------------------- ----------------------
3 4 ---------------------- ----------------------
3 8 ---------------------- ----------------------
3 16 ---------------------- ----------------------
In quadrature, we know we're only approximating an integral. It is often useful to have an estimate of the error made in approximation. While this estimate can't be taken completely seriously, it is often used as a guide to deciding when an approximation is probably accurate enough, or when it is likely that using more points would produce a better answer.
To approximate the error, we simply make two estimates of the integral, in such a way that we expect the second estimate to be significantly more accurate. We could, for instance compare two rules of different order, or two composite rules, one of which uses twice as many subintervals. In both cases, the difference between the two approximations to the integral is an estimate of the error.
Exercise: consider the Runge function 1/(1+x^2) on the interval [-4,4]. The exact integral is 2*atan(4). Compute the MATLAB Sum estimate for the integral, using the following number of points. Use the difference of the first and second sums as an error estimate for the first sum. Use the difference of the second and third sums to estimate the error in the second sum, and so on.
N MATLAB Sum Estimated Error True Error
2 ----------------- -------------- -----------
5 ----------------- -------------- -----------
17 ----------------- -------------- -----------
65 ----------------- -------------- -----------
257 -----------------
Using any method discussed in this lab, estimate the integral over [-1,1] of
0.92 * cosh ( x ) - cos ( x )Mail me the name of the quadrature method you used, the order N, the number of subintervals M, and the estimated value of the integral. Your integral estimate must be correct in the first three significant digits. Please turn your results in by Friday, December 10.