|Matlab hint||Exercise 2|
|The Midpoint Method||Exercise 3|
|The Trapezoid Method||Exercise 5|
|Singular Integrals||Exercise 6|
|Newton-Cotes Rules||Exercise 7|
|Gauss Quadrature||Exercise 8|
|Adaptive quadrature||Exercise 9|
The term ``numerical quadrature'' refers to the estimation of an area, or, more generally, any integral. (You might hear the term ``cubature'' referring to estimating a volume in three dimensions, but most people use the term ``quadrature.'') We might want to integrate some function 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 first discuss the ``degree of precision'' of a quadrature rule, and its relation to the ``order of accuracy.'' We consider some simple tests to measure the degree of precision and the order of accuracy 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 either increasing the order of the rule or subdividing the interval. Finally, we will consider a way of adapting the details of the integration method to meet a desired error estimate. In the majority of the exercises below, we will be more interested in the error (difference between calculated and known value) than in the calculated value itself.
A word of caution. We discuss three similar-sounding concepts:
This lab will take four sessions. If you print this lab, you may prefer to use the pdf version.
Matlab provides the capability of defining functions ``in line'' instead
of writing m-files to do it. This feature is most convenient when the
function to be defined is very simple-a line of code, say-or when
you have a function that requires several arguments but you only
are interested in varying one of them. The way to do this is to
use the inline function. (The full help description is
available from the command line using help inline, from
the help facility on your PC, or from
the Mathworks on-line reference.)
Suppose, for example, you want to define a function
You could do this by writing
sq=inline('x.^2')(Letters such as x are understood to mean independent variables.) or, less ambiguously, by writing
sq=inline('x.^2','x')You could then use sq(x) later, just as if you had defined it in an m-file. Alternatively, you could use it where a function name is needed. In the next section you will be writing an integration function named midpoint that requires a function name as its first argument. If you wanted to apply it to the integral , you might write
There is a nice way to use inline to streamline a sequence of calculations computing the integrals of ever higher degree polynomials in order to find the degree of precision of a quadrature rule. The following statement
q=midpointquad(inline('5*x.^4','x'),0,1,11);1-qfirst computes using the midpoint rule, and then prints the error (=1-q because the exact answer is 1). You would only have to change
'4*x.^3'to check the error in , and you can make that change with judicious use of the arrow keys.
A second remark is that errors should be reported in scientific notation (like 1.234e-3, not .0012). This is particularly important if you want to visually estimate ratios of errors. In this lab, I would much prefer that you use full-precision values when you compute ratios of errros. For example, you might use code like
err20=midpointquad('runge',-5,5,20)-2*atan(5); err40=midpointquad('runge',-5,5,40)-2*atan(5); ratio=err20/err40will yield a ratio of errors without loss of accuracy due to reading numbers off the computer screen.
In general, numerical quadrature involves breaking an interval into subintervals, estimating or modelling the function on each subinterval and integrating it there, then adding up the partial integrals.
Perhaps the simplest method of numerical integration is the midpoint method (Atkinson, p. 269). This method is based on interpolation of the integrand by the constant and multiplying by the width of the interval. The result is a form of Riemann sum that you probably learned in elementary calculus when you first studied integration.
Break the interval into subintervals with endpoints (there is one more endpoint than intervals, of course). Then the midpoint rule can be written as
function quad = midpointquad( f, a, b, n) % comments % your name and the datewhere f indicates the name of a function, a and b are the lower and upper limits of integration, and n is the number of points, not the number of intervals. The code for your m-file might look like the following:
xpts = linspace( ??? ) ; h = ??? ; % length of subintervals xmidpts = 0.5 * ( xpts(1:n-1) + xpts(2:n) ); fmidpts = ??? quad = h * sum ( fmidpts );
n h Midpoint Result Error 11 1.0 _________________ __________________ 101 0.1 _________________ __________________ 1001 0.01 _________________ __________________ 10001 0.001 _________________ __________________
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. (Atkinson, page 266.)
To determine the degree of precision of a rule, we might look at
the approximations of the integrals
f Midpoint Result Error 1 ___________________ ___________________ 2 * x ___________________ ___________________ 3 * x^2 ___________________ ___________________ 4 * x^3 ___________________ ___________________
The trapezoid rule breaks [a,b] into subintervals, approximates the integral on each subinterval as the product of its width times the average function value, and then adds up all the subinterval results, much like the midpoint rule. The difference is in how the function is approximated. The trapezoid rule can be written as
To apply the trapezoid rule, we need to generate points and evaluate the function at each of them. Then, apply either (2) or (3) as appropriate.
function quad = trapezoidquad( f, a, b, n ) % comments % your name and the dateYou may use either form of the trapezoid rule.
f Trapezoid Result Error 1 ___________________ ___________________ 2 * x ___________________ ___________________ 3 * x^2 ___________________ ___________________ 4 * x^3 ___________________ ___________________
n h Trapezoid Result Error 11 1.0 _________________ __________________ 101 0.1 _________________ __________________ 1001 0.01 _________________ __________________ 10001 0.001 _________________ __________________
The midpoint and trapezoid rules seem to have the same precision and about the same accuracy. There is a difference between them, though. Some integrals have perfectly well-defined values even though the integrand has some sort of mild singularity. Atkinson discusses some sophisticated ways to perform these integrals, but there is a simple way that can be adequate for the case that the singularity appears at the endpoint of an interval. Something is lost, however.
Consider the integral
n h Midpoint Result Error 11 0.1 _________________ __________________ 101 0.01 _________________ __________________ 1001 0.001 _________________ __________________ 10001 0.0001 _________________ __________________Estimate the rate of convergence (power of ) as . You should see that the singularity causes a loss in the rate of convergence.
Look at the trapezoid rule for a minute. One way of interpreting that rule is to say that if the function is roughly linear over the subinterval , then the integral of is the integral of the linear function that agrees with (i.e., interpolates ) over that interval. What about trying higher order methods? It turns out that Simpson's rule can be derived by picking triples of points, interpolating the integrand by a quadratic polynomial, and integrating the quadratic. The trapezoid rule and Simpson's rule are Newton-Cotes rules of index one and index two, respectively. In general, a Newton-Cotes formula uses the idea that if you approximate a function by a polynomial interpolant, then you can approximate the integral of that function with the integral of the polynomial interpolant. The polynomial interpolant in this case being taken on a uniformly distributed set of points, including the end points. Atkinson defines the ``index'' of a Newton-Cotes rule as one fewer than the number of points it uses.
We applied the trapezoid rule to an interval by breaking it into subintervals and repeatedly applying a simple formula for the integral on a single subinterval. Similarly, we will be constructing higher-order rules by repeatedly applying Newton-Cotes rules over subintervals. But Newton-Cotes formulæ are not so simple as the trapezoid rule, so we will first write a helper function to apply the rule on a single subinterval.
Over a single interval, all (closed) Newton-Cotes formulæ can be written as
function quad = nc_single ( f, a, b, index ) % comments % your name and the dateThere are no subintervals in this case. The coding might look like something like this:
xvec = linspace ( a, b, index+1 ); wvec = nc_weight ( index ); fvec = ??? quad = (b-a) * dot(wvec , fvec);
f Error Error Error index=3 index=4 index=5 4 * x^3 __________ __________ ___________ 5 * x^4 __________ __________ ___________ 6 * x^5 __________ __________ ___________ 7 * x^6 __________ __________ ___________ Degree ___ ___ ___
The objective of numerical quadrature rules is to accurately approximate integrals. We have already seen that polynomial interpolation on uniformly spaced points does not always converge, so it should be no surprise that increasing the order of Newton-Cotes integration might not produce accurate quadratures.
index nc_single Result Error 3 _________________ __________________ 7 _________________ __________________ 11 _________________ __________________ 15 _________________ __________________
The results of the previous exercise should have convinced you that you need a composite Newton-Cotes rule in order to accurately do numerical integration. In the following exercise you will use nc_single as a helper function for a composite Newton-Cotes routine.
function quad = nc_quad( f, a, b, index, numSubintervals) % commentsThis function will perform these steps: (1) break the interval into numSubintervals subintervals; (2) use nc_single to integrate over each subinterval; and, (3) add them up.
nc_quad('runge', -5, 5, 3, 10) = 2.74533025
Subin- tervals index nc_quad Error Err ratio 10 1 ____________ _____________ __________ 20 1 ____________ _____________ __________ 40 1 ____________ _____________ __________ 80 1 ____________ _____________ __________ 160 1 ____________ _____________ __________ 320 1 ____________ _____________ 10 2 ____________ _____________ __________ 20 2 ____________ _____________ __________ 40 2 ____________ _____________ __________ 80 2 ____________ _____________ __________ 160 2 ____________ _____________ __________ 320 2 ____________ _____________
In the previous exercise, the table served to illustrate the behavior of the integration routine. Suppose, on the other hand, that you had an integration routine and you wanted to be sure it had no errors. It is not good enough to just see that you can get ``good'' answers. In addition, it must converge at the correct rate. I use tables such as the previous one to convince myself that newly-written code is correct.
Like Newton-Cotes quadrature, Gauss-Legendre quadrature interpolates the integrand by a polynomial and integrates the polynomial. Instead of uniformly spaced points, Gauss-Legendre uses optimally-spaced points. Furthermore, Gauss-Legendre converges as degree gets large, unlike Newton-Cotes, as we saw above. Of course, in real applications, one does not use higher and higher degrees of quadrature-one uses more and more subintervals, each with some given degree of quadrature.
The disadvantage of Gauss-Legendre quadrature is that there is no simple formula for the node points and weights. Instead, the values are tabulated (see, for example, Table 5.10 in Atkinson, page 276). We will be using a Matlab function to serve as a table of node points and weights. Atkinson discusses Gauss-Legendre quadrature in Section 5.3.
Normally, Gauss-Legendre quadrature is characterized by the number of integration points. We speak of three-point Gauss, etc. In this lab, we will regard this as the ``index'' in order to be consistent with Newton-Cotes quadrature and also to distinguish the number of integration points from the number of subintervals in composite Gauss-Legendre quadrature.
The following two exercises involve writing m-files analogous to nc_single.m and nc_quad.m.
function quad = gl_single ( f, a, b, index ) % commentsAs with nc_single there are no subintervals in this case. Your coding might look like something like this:
[xvec, wvec] = gl_weight ( a, b, index ); fvec = ??? quad = dot( wvec , fvec );
f Error Error index=2 index=3 3 * x^2 __________ ___________ 4 * x^3 __________ ___________ 5 * x^4 __________ ___________ 6 * x^5 __________ ___________ 7 * x^6 __________ ___________ Degree ___ ___
index gl_single Result Error 3 _________________ __________________ 7 _________________ __________________ 11 _________________ __________________ 15 _________________ __________________
You might be surprised at how much better Gauss-Legendre integration is than Newton-Cotes, using a single interval. There is a similar advantage for composite integration, but it is hard to see for small indices. 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 is compounded 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.
function quad = gl_quad( f, a, b, index, numSubintervals) % commentsThis function will perform two steps: (1) break the interval into numSubintervals subintervals; (2) use gl_single to integrate over each subinterval; and, (3) add them up.
gl_quad('runge', -5, 5, 4, 10) = 2.7468113
Subin- tervals index gl_quad Error Err ratio 10 1 ____________ _____________ __________ 20 1 ____________ _____________ __________ 40 1 ____________ _____________ __________ 80 1 ____________ _____________ __________ 160 1 ____________ _____________ __________ 320 1 ____________ _____________ 10 2 ____________ _____________ __________ 20 2 ____________ _____________ __________ 40 2 ____________ _____________ __________ 80 2 ____________ _____________ __________ 160 2 ____________ _____________ __________ 320 2 ____________ _____________
Our final task will consider ``adaptive quadrature.'' Adaptive quadrature employs non-uniform division of the interval of integration into subintervals of non-equal length. It uses smaller subintervals where the integrand is changing rapidly and larger subintervals where it is flatter. The advantage of this approach is that it minimizes the work necessary to compute a given integral. Atkinson discusses adaptive quadrature starting on page 300. There, he expresses it as a recursive algorithm and we will follow his lead here. (This is a case where the potential inefficiency of recursive algorithms is greatly overshadowed by their simplicity. You might want to think how you might recast the recursive function below using loops.)
Numerical integration is used often with integrands that are very complicated and take a long time to compute. Although this section will use simple integrands for illustrative purposes, you should think of each evaluation of the integrand (``function call'') to take a long time. Thus, the objective is to reach a given accuracy with a minimum number of function calls.
Recall that, if an integration method has degree of precision , then the local error on an integration interval of length can be estimated by dividing the interval into two subintervals of length each. Denote by the integral over the interval of length and and the two estimates of the integral on the left and right subintervals. Then
The basic structure of one simple adaptive algorithm depends on using the error estimate (6) over each integration subinterval. If the error is acceptably small over each interval, the process stops, and, if not, continues recursively. In the following exercise, you will write a recursive function to implement this procedure.
function [Q,errEst,x,recursions]= ... adaptquad(f,x0,x1,tol,recursions) % [Q,errEst,x,recursions]= % adaptquad(f,x0,x1,tol,recursions) % adaptive quadrature % input parameters % f = function to integrate % x0 = left end point % x1 = right end point % tol = desired accuracy % recursions = number of allowable recursions left % % output parameters % Q = estimate of the value of the integral % errEst = estimate of error in Q % x = all intermediate integration points % recursions = minimum number of recursions % remaining after convergence % Add a mid-point and re-estimate integral xmid=(x0+x1)/2; % Qleft and Qright are integrals over two halves INDEX=3; Qboth=gl_single(f,x0,x1,INDEX); Qleft=gl_single(f,x0,xmid,INDEX); Qright=gl_single( ??? ); % p=degree of precision of Gauss-Legendre p=2*INDEX-1; errEst= ??? ; if errEst<tol | recursions<=0 %vertical bar means "or" % either ran out of recursions or converged Q= ??? ; x=[x0 xmid x1]; else % not converged -- do it again [Qleft,estLeft,xleft,recursLeft]=adaptquad(f, ... x0,xmid,tol/2,recursions-1); [Qright,estRight,xright,recursRight]=adaptquad(f, ... ??? ); % recursive work is all done, return answers % don't want xmid to appear twice in x x=[xleft xright(2:length(xright))]; Q= ??? ; errEst= ??? ; recursions=min(recursLeft,recursRight); endNote: The input and output parameter recursions is not theoretically necessary, but is used to guard against infinite recursion. The output vector x is not necessary, either, but will be used to show the effect of the adaptation.
adaptquad for Runge tol est. error exact error 1.e-3 __________ __________ 1.e-6 __________ __________ 1.e-9 __________ __________You should find that the estimated and exact errors are close in size, and smaller than tol.
Now you are ready to look a little more closely at the results of this recursive adaptive algorithm. Some of the questions you might ask are listed below.
In the following exercises, you will examine two functions that are more difficult to integrate. The first is a scaled version of the Runge function, , where The value of the integral on [-1,1] of the scaled Runge function is
The second function is , and the value of its integral over the interval [-1,1] is
A third function that is difficult to integrate is . This function has an integrable singularity at that is close to being nonintegrable.
[Q,estErr,x,recursions]=adaptquad('srunge',-1,1,1.e-10,50);What are the estimated and true errors?
xave=(x(2:length(x))+x(1:length(x)-1))/2; dx= x(2:length(x))-x(1:length(x)-1); semilogy(xave,dx,'*')A semilog plot is appropriate here because of the wide range of interval sizes. Please include this plot with your summary.
xave=(x(2:length(x))+x(1:length(x)-1))/2; dx= x(2:length(x))-x(1:length(x)-1); semilogy(xave,dx,'*')where x is replaced by the variable name you used. A semilog plot is appropriate here because of the wide range of interval sizes. Please include this plot with your summary.
In the following exercise you will apply adaptquad to the function with an integral singularity in Exercise 4. In this example, you will see the benefit of the recursions variable. This variable takes the value 0 when the convergence test has failed too many times.
Back to MATH2070 page.