|
|
||
|
|
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. Then 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.
In general, numerical quadrature involves breaking an interval [a,b] into subintervals, sometimes of equal size and sometimes not. The midpoint rule breaks [a,b] into equal subintervals, approximates the integral one each subinterval as the product of its width h times the function value at the midpoint, and then adds up all the subinterval results.
Midpoint Method
Divide the interval [a.b] into (n-1) subintervals [xi,xi+1]
Repeat steps 3 through 5 for each subinterval
hi=(xi+1-xi)
xi+1/2=(xi+xi+1)/2
Qi=hi*f(xi+1/2)
Q=sum of Qi
It's a little tricky to do a composite midpoint rule with MATLAB, because we need to compute the locations of the midpoints xi+1/2. The easiest way to compute the midpoints is to compute the n endpoints into a vector of points and then average pairs of them, giving a new vector of n-1 midpoints. The vector of midpoints will contain xi+1/2 in position i.
Write a routine called mid_quad.m, which is declared as:
function Q = mid_quad ( f, a, b, n )
The coding might look like:
h = ( b - a ) / ( n - 1 ) % all subintervals are the same size
xpts = ? % n evenly spaced points, including endpoints
xmids = 0.5 * ( xpts(1:n-1) + xpts(2:n) ) % midpoints of intervals. length(xmids) is n-1
fvec = ?
Q = h * sum ( fvec )
Exercise 1: Use mid_quad to estimate the integral of 5*x4 over [0,1], using the given values of n, and record the error, noting that the correct answer is 1.0.
n h Midpoint Result Error
2 1 ---------------------- ----------------------
11 0.1 ---------------------- ----------------------
101 0.01 ---------------------- ----------------------
1001 0.001 ---------------------- ----------------------
How fast is the error decreasing with h?
Matlab allows you to string several commands on a single line, separated by semicolons. This feature makes it easy to do an exercise such as the Exercise 1 with the following statement:
n=2;q=mid_quad(inline('5*x.^4'),0,1,n);[n,q,1-q]
and then use the uparrow key to recall the statement, change n, and do it again.
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) = 1
p1(x) = 2 * x
p2(x) = 3 * x2
...
pk(x) = ( k + 1 ) * xk
since the exact value of every one of these integrals is 1.
Exercise 2: 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 exact 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 trapezoid rule breaks [a,b] into subintervals, approximates the integral one each subinterval as the product of its width h times the average function value, and then adds up all the subinterval results.
Trapezoid Method
Divide the interval [a.b] into (n-1) subintervals [xi,xi+1]
Repeat steps 3 and 4 for each subinterval
hi=(xi+1-xi)
Qi=hi*( f(xi) + f(xi+1) )/2
Q=sum of Qi
To apply the trapezoid rule, we need to generate N points and evaluate the function at each of them. Then these values are used to evaluate Q using the above formula for each interval, and then add them up. Note that, except for the very ends, the (1/2)'s add up to 1's. Write out he sum for N=4 on a piece of paper to see what is happening.
Write a routine called trap_quad.m, that computes the trapezoid rule using uniform mesh increments and is declared as:
function Q = trap_quad ( f, a, b, n )
The coding might look like something like this:
h = ( b - a ) / ( n - 1 )
xvec = ?
fvec = ?
Q = h * ( sum ( fvec ) - 0.5 * ( fvec(1) + fvec(n) ) )
Exercise 3: 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?
Exercise 4: To study the precision of the trapezoid rule, use a single interval ( n = 2 ), and estimate the integrals of the same test functions used for the rectangle rule over [0,1]. The exact 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 trapezoid rule?
Look at the trapezoid rule for a minute. One way of interpreting that rule is to say that if the function f is roughly linear over the subinterval [xi,xi+1], then the integral of f is the integral of the linear function that agrees with f (i.e., interpolates f) over that interval. What about trying higher order methods? It turns out that Simpson's rule is derived by picking triples of points, interpolating the integrand f by a quadratic polynomial, and integrating the quadratic. These two rules are Newton-Cotes rules of order two and three, respectively. In general, 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.
Write a routine called nc_single.m, which is declared as:
function Q = nc_single ( f, a, b, n )
Be careful here because the n refers to the order. There are no subintervals in this case. The coding might look like something like this:
h = ( b - a ); %This is not divided by (n-1) because we are doing a single interval.
xvec = linspace ( a, b, n );
wvec = nc_weight ( n );
fvec = ?
Q = h * dot(wvec , fvec);
The routine nc_weight.m returns the Newton-Cotes coefficients for the rule of given order. It is available from the web page.
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. 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 5: Consider the Runge function 1/(1+x^2) on the interval [-5,5]. The exact integral is 2*atan(5). Use rules of increasing order, each time recording the quadrature error.
N Newton-Cotes Result Error
3 ---------------------- ----------------------
5 ---------------------- ----------------------
9 ---------------------- ----------------------
17 ---------------------- ----------------------
Using more precise Newton-Cotes rules means using increasing degree polynomial interpolation with equally spaced nodes. This wasn't successful in earlier labs, it should be no surprise that it doesn't work here, either.
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, nc_comp.m, which looks like
function Q = nc_comp ( f, a, b, norder , nsubints )
which approximates the integral of f(x) over [A,B] using a Newton-Cotes rule of order norder, breaking the interval into nsubints subintervals. The routine might look like this:
endpts = linspace ( a, b, nsubints+1)
Q = 0
for i = 1:nsubints
Q = Q + nc_single ( f, endpts(i), endpts(i+1), norder )
Note that we are computing nsubints+1 subinterval endpoints, which define nsubints subintervals. Look carefully at the code: there are nsubints intervals, each of which is divided into norder-1 intervals, for a total of nsubints*(norder-1) intervals.
Exercise 6: Consider the Runge function 1/(1+x^2) on the interval [-5,5]. The exact integral is 2*atan(5). Try the composite rules of order norder, with nsubints intervals, each time recording the error.
norder nsubints Composite Newton-Cotes Error Total number of intervals
2 4 ---------------------- --------------- -------
2 8 ---------------------- --------------- -------
2 16 ---------------------- --------------- -------
4 4 ---------------------- --------------- -------
4 8 ---------------------- --------------- -------
4 16 ---------------------- --------------- -------
8 2 ---------------------- --------------- -------
8 4 ---------------------- --------------- -------
8 8 ---------------------- --------------- -------
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.
The down side of Gauss's quadrature methods (more properly called Gauss-Legendre methods) is that there is no simple formula for the weights and node points. Instead, people used to have tables of these values in books when they did their work by hand. Today, computer programs are available for computing these quantities. For our purposes, the weights will be computed by the routine gl_weight.m, and the location of the points by gl_space.m, both available from the web page.
(Gauss's points are actually the places where the Legendre polynomial is zero. It is not required, but you may wish to 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 gl_quad.m, which is declared as:
function Q = gl_quad ( 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 = ?
Q = h * dot( wvec , fvec )
Exercise 7: Consider the Runge function 1/(1+x^2) on the interval [-5,5]. The exact integral is 2*atan(5). Use rules of increasing order, each time recording the quadrature error. You can tell how many x points are used for evaluation by looking at gl_space.m. See for yourself how much more efficient Gauss-Legendre quadrature is than Newton-Cotes.
N Gauss-Legendre Result Error Number of evaluation points
2 ---------------------- ----------- -----------
4 ---------------------- ----------- -----------
8 ---------------------- ----------- -----------
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. In fact, when Gauss-Legendre integration is used in a computer program, it is generally in the form of a composite formulation because it is difficult to compute the weights and integration points accurately for high order Gauss-Legendre integration. The efficiency of Gauss-Legendre integration becomes even greater in multiple dimensions, and essentially all computer programs that use the finite element method use composite Gauss-Legendre integration rules to compute the coefficient matrices.
In quadrature, we're approximating an integral--a scalar value. It is often useful to have an estimate of the error made in approximation. While this estimate can't be taken completely seriously (there is no proof), 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 8: Consider the Runge function 1/(1+x^2) on the interval [-5,5]. The exact integral is 2*atan(5). Estimate the integral using the mid_quad.m file (midpoint rule) you wrote for Exercise 1, with the following numbers of points. Use the difference of the first and second estimates 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 Estimated value Estimated Error True Error
2 ----------------- -------------- -----------
5 ----------------- -------------- -----------
17 ----------------- -------------- -----------
65 ----------------- -------------- -----------
257 ----------------- -----------
Exercise 9: Using any quadrature method discussed in this lab, estimate to three significant digits the integral over [-1,1] of
0.92 * cosh ( x ) - cos ( x )
Why do you think your estimate is correct to three significant digits?
Back to the MATH2070 page.
Last revised on November 26, 2000.