Unformatted Document Excerpt
Coursehero >>
Pennsylvania >>
Pittsburgh >>
MATH 2070
Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.
Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.
LAB MATH2070: 9: Quadrature
Introduction Matlab hint The Midpoint Method Precision The Trapezoid Method Singular Integrals Newton-Cotes Rules Gauss Quadrature Adaptive quadrature
Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise Exercise
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1
Introduction
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 f (x) or a set of tabulated data. The domain might be a nite or innite interval, it might be a rectangle or irregular shape, it might be a three dimensional cell. We rst 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 midpointquad and trapezoidal rules. Then we consider the NewtonCotes 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 (dierence between calculated and known value) than in the calculated value itself. A word of caution. We discuss three similar-sounding concepts: Degree of precision: the largest value of n so that all polynomials of degree n and below are integrated exactly. (Degree of a polynomial is the highest power of x appearing in it.) Order of accuracy: the value of n so that the error is O(hn ), where h measures the subinterval size. Index: a number distinguishing one of a collection of rules from another. These can be related to one another, but are not the same thing. This lab will take four sessions. If you print this lab, you may prefer to use the pdf version.
2
Matlab hint
Matlab provides the capability of dening functions in line instead of writing m-les to do it. This feature is most convenient when the function to be dened is very simplea line of code, sayor when you have a 1
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 dene a function sq(x)=x^2. 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 dened it in an m-le. 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 rst argument. If you wanted to apply it to the integral 1 2 x dx, you might write 0 q=midpointquad(inline(x.^2,x),0,1,11) 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 nd the degree of precision of a quadrature rule. The following statement q=midpointquad(inline(5*x.^4,x),0,1,11);1-q rst computes 0 5x4 dx using the midpoint rule, and then prints the error (=1-q because the exact answer 1 is 1). You would only have to change 5*x.^4 into 4*x.^3 to check the error in 0 4x3 dx, and you can make that change with judicious use of the arrow keys. A second remark is that errors should be reported in scientic notation (like 1.234e-3, not .0012). You can force Matlab to display numbers in this format using the command format short e (or format long e for 15 decimal places). This is particularly important if you want to visually estimate ratios of errors. In this lab, you should 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/err40 to get a ratio of errors without loss of accuracy due to reading numbers o the computer screen.
1
3
The Midpoint Method
In general, numerical quadrature involves breaking an interval [a, b] 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 (Ralston and Rabinowitz, p. 120, or Atkinson, p. 269). This method is based on interpolation of the integrand f (x) by the constant f ( a+b ) and multiplying by the width of the interval. The result is a form of Riemann sum that you probably 2 saw in elementary calculus when you rst studied integration. Break the interval [a, b] into N 1 subintervals with endpoints x1 , x2 , . . . , xN 1 , xN (there is one more endpoint than intervals, of course). Then the midpoint rule can be written as
N 1
Midpoint rule =
k=1
(xk+1 xk )f (
xk + xk+1 ). 2
(1)
In the exercise that follows, you will be writing a Matlab function to implement the midpoint rule. 2
Exercise 1: (a) Write a function m-le called midpointquad.m with signature function quad = midpointquad( f, a, b, n) % quad = midpointquad( f, a, b, n) % comments % your name and the date where 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-le 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 ); (b) Test your midpointquad routine by computing 0 2xdx = 1. Even if you use only one interval (i.e. n=2) you should get the exact answer because the midpoint rule integrates linear functions exactly. (c) Use your midpoint routine to estimate the integral of our friend, the Runge function, f (x) = 1/(1 + x2 ), over the interval [5, 5]. (If you do not have a copy of the Runge function handy, you can download my version of runge.m.) The exact answer is 2*atan(5). Fill in the following table, using scientic notation for the error values so you can see the pattern. n 11 101 1001 10001 h 1.0 0.1 0.01 0.001 Midpoint Result _________________ _________________ _________________ _________________ Error __________________ __________________ __________________ __________________
1
(d) Estimate the order of accuracy (an integer power of h) by examining the behavior of the error when h is divided by 10. (In previous labs, we have estimated such orders by repeatedly doubling the number of subintervals. Here, we multiply by ten. The idea is the same.)
4
Precision
If a quadrature rule can compute exactly the integral of any polynomial up to some specic degree, we will call this its degree of precision. Thus a rule that can correctly integrate any cubic, but not quartics, has precision 3. (This term is not explicitly used by Ralston and Rabinowitz, but the concept is used. The term is used by Atkinson, page 266.)
3
To determine the degree of precision of a rule, we might look at the approximations of the integrals
1
1dx
0 1
= = = . . .
[x]1 0 [x2 ]1 0 [x3 ]1 0
2xdx
0 1 0
3x2 dx
1 0
(k + 1)xk dx
=
[xk+1 ]1 0
The exact value of each one of these integrals is 1. Exercise 2: (a) To study the degree of precision of the midpoint method, use a single interval (i.e. n = 2), and estimate the integrals of the test functions over [0,1]. The exact answer is 1 each time. f Midpoint Result Error ___________________ ___________________ ___________________ ___________________
1 ___________________ 2 * x ___________________ 3 * x^2 ___________________ 4 * x^3 ___________________
(b) What is the degree of precision of the midpoint rule? (c) For some methods, but not all, the degree of precision is one less than the order of accuracy. Is that the case for the midpoint rule?
5
The Trapezoid Method
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 dierence is in how the function is approximated. The trapezoid rule can be written as
N 1
Trapezoid rule =
k=1
(xk+1 xk )
f (xk ) + f (xk+1 ) 2
(2)
If you compare the midpoint rule (1) and the trapezoid rule, you will see that the midpoint rule takes f at the midpoint of the subinterval and the trapezoid takes the average of f at the endpoints. If each of the subintervals happens to have length h, then the trapezoid rule becomes h h f (x1 ) + f (xN ) + h 2 2
N 1
f (xk ).
k=2
(3)
To apply the trapezoid rule, we need to generate N points and evaluate the function at each of them. Then, apply either (2) or (3) as appropriate. Exercise 3: 4
(a) Use your midpointquad.m m-le as a model and write a function m-le called trapezoidquad.m to evaluate the trapezoid rule. The signature of your m-le should be function quad = trapezoidquad( f, a, b, n ) % quad = trapezoidquad( f, a, b, n ) % more comments % your name and the date You may use either form of the trapezoid rule. (b) To test your routine and 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 midpoint rule over [0,1]. The exact answer should be 1 each time. f Trapezoid Result Error ___________________ ___________________ ___________________ ___________________
1 ___________________ 2 * x ___________________ 3 * x^2 ___________________ 4 * x^3 ___________________
(c) What is the degree of precision of the trapezoid rule? (d) Use the trapezoid method to estimate the integral of the Runge function over [5, 5], using the given values of n, and record the error using scientic notation. n 11 101 1001 10001 h 1.0 0.1 0.01 0.001 Trapezoid Result _________________ _________________ _________________ _________________ Error __________________ __________________ __________________ __________________
(e) Estimate the rate at which the error decreases as h decreases. (Find the power of h that best ts the error behavior.) This is the order of accuracy of the rule. (f) For some methods, but not all, the degree of precision is one less than the order of accuracy. Is that the case for the trapezoid rule?
6
Singular Integrals
The midpoint and trapezoid rules seem to have the same precision and about the same accuracy. There is a dierence between them, though. Some integrals have perfectly well-dened values even though the integrand has some sort of mild singularity. There are 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 1 log(x)dx = /2, I=
0
where log refers to the natural logarithm. Note that the integrand is innite at the left endpoint, so you could not use the trapezoid rule to evaluate it. The midpoint rule, conveniently, does not need the endpoint values. Exercise 4: Apply the midpoint rule to the above integral, and ll in the following table.
5
n 11 101 1001 10001
h 0.1 0.01 0.001 0.0001
Midpoint Result _________________ _________________ _________________ _________________
Error __________________ __________________ __________________ __________________
Estimate the rate of convergence (power of h) as h 0. You should see that the singularity causes a loss in the rate of convergence.
7
Newton-Cotes Rules
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 [xk , xk+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 Simpsons rule can be derived by picking triples of points, interpolating the integrand f by a quadratic polynomial, and integrating the quadratic. The trapezoid rule and Simpsons 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. This idea does not work very well for derivatives. The polynomial interpolant in this case being taken on a uniformly distributed set of points, including the end points. Ralston and Rabinowitz do not use the term index, but the variable n that is used to characterize dierent formul corresponds to the index and is one fewer than the number of points in the formula. Atkinson explicitly denes 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 rst 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
b n
f (x)dx Qn (f ) =
a k=0
wk,n f (xk )
where f is a function and xk are n + 1 evenly-spaced points between a and b. The weights wk,n can be computed from the Lagrange interpolation polynomials k,n as
1
wk,n = (b a)
0
k,n ()d.
(The Lagrange interpolation polynomials arise because we are doing a polynomial interpolation. See Ralston and Rabinowitz, p. 118 or Atkinson, p. 263.) The weights do not depend on f , and on a and b in a simple manner, so they are often tabulated for the unit interval. In the exercise below, I will provide them through a function. Exercise 5: (a) Download nc weight.m. (b) Write a routine called nc single.m with the signature function quad = nc_single ( f, a, b, index ) % quad = nc_single ( f, a, b, index ) 6
% more comments % your name and the date There are no subintervals in this case. The coding might look like something like this: xvec wvec fvec quad = = = = linspace ( a, b, index+1 ); nc_weight ( index ); ??? (b-a) * sum(wvec .* fvec);
1 0
(c) Test your function by showing its precision is at least 1 for index 1:
2xdx = 1 exactly.
(d) Fill in the following table by computing the integrals over [0,1] of the indicated integrands using nc single. Ralston and Rabinowitz (p. 119) and Atkinson (p. 266) indicate that the degree of precision is equal to the index when the index is odd and the degree of precision is one more than the index when the index is even. Your results should agree, further conrming that your function is correct. (Hint: You can use inline to simplify your work.) f 4 * x^3 5 * x^4 6 * x^5 7 * x^6 Degree Error index=3 __________ __________ __________ __________ ___ Error index=4 __________ __________ __________ __________ ___ Error index=5 ___________ ___________ ___________ ___________ ___
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. Exercise 6: Attempt to get accurate estimates of the integral of the Runge function over the interval [-5,5]. Recall that the exact answer is 2*atan(5). Fill in the following table index 3 7 11 15 nc_single Result _________________ _________________ _________________ _________________ Error __________________ __________________ __________________ __________________
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. Exercise 7: (a) Write a function m-le called nc quad.m to perform a composite Newton-Cotes integration. Use the following signature. function quad = nc_quad( f, a, b, index, numSubintervals) % quad = nc_quad( f, a, b, index, numSubintervals) % comments % your name and the date 7
This 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. (b) The most elementary test to make when you write this kind of routine is to check that you get the same answer when numSubintervals=1 as you would have obtained using nc single. Choose at least one line from the table in Exercise 5 and make sure you get the same result using nc quad. (c) Test your routine by checking the following value nc_quad(runge, -5, 5, 3, 10) = 2.74533025 (d) Fill in the following table using the Runge function. Subintervals index 10 20 40 80 160 320 10 20 40 80 160 320 10 20 40 80 160 320 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 nc_quad Error _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ _____________ Err ratio __________ __________ __________ __________ __________
__________ __________ __________ __________ __________
__________ __________ __________ __________ __________ __________
(e) For each index, estimate the order of convergence by taking the sequence of ratios of the error for N subintervals divided by the error for (2*N) subintervals and guessing the power of two that best approximates the limit of the sequence. the In 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. Tables such as the previous one are one of the most powerful debugging and verication tools a researcher has.
8
Gauss Quadrature
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 quadratureone uses more and more subintervals, each with some given degree of quadrature. 8
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, Ralston and Rabinowitz, Table 4.1, p. 100 or 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 NewtonCotes 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-les analogous to nc single.m and nc quad.m. Exercise 8: (a) Download the le gl weight.m. This le returns both the node points and weights for GaussLegendre quadrature for index=n points. (b) Write a routine called gl single.m with the signature function quad = gl_single ( f, a, b, index ) % quad = gl_single ( f, a, b, index ) % comments % your name and the date As 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 = sum( wvec .* fvec ); (c) Test your function by showing its precision is at least 1 for one interval:
1 0
2xdx = 1 exactly.
(d) Fill in the following table by computing the integrals over [0,1] of the indicated integrands using gl single. Ralston and Rabinowitz (page 98-99) and Atkinson (page 272) show that the degree of precision of the method is 2n 1, where n means index, and your results should agree, further conrming that your function is correct. (Hint: You can use inline to simplify your work.) f 3 * x^2 4 * x^3 5 * x^4 6 * x^5 7 * x^6 Degree Error index=2 __________ __________ __________ __________ __________ ___ Error index=3 ___________ ___________ ___________ ___________ ___________ ___
(e) Get accurate estimates of the integral of the Runge function over the interval [-5,5]. Recall that the exact answer is 2*atan(5). Fill in the following table index 3 7 11 15 gl_single Result _________________ _________________ _________________ _________________ Error __________________ __________________ __________________ __________________
9
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 dicult to compute the weights and integration points accurately for high order Gauss-Legendre integration. The eciency of Gauss-Legendre integration is compounded in multiple dimensions, and essentially all computer programs that use the nite element method use composite Gauss-Legendre integration rules to compute the coecient matrices. Exercise 9: (a) Write a function m-le called gl quad.m to perform a composite Gauss-Legendre integration. Use the following signature. function quad = gl_quad( f, a, b, index, numSubintervals) % quad = gl_quad( f, a, b, index, numSubintervals) % comments % your name and the date This 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. (b) The most elementary test to make when you write this kind of routine is to check that you get the same answer when numSubintervals=1 as you would have obtained using gl single. Choose at least one line from the table in the previous exercise (8) and make sure you get the same result using gl quad. (c) Test your routine by checking the following value gl_quad(runge, -5, 5, 4, 10) = 2.7468113 (d) Fill in the following table. Subingl_quad tervals index Error 10 20 40 80 160 320 10 20 40 80 160 320 45 90 46 92 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________ Err ratio __________ __________ __________ __________ __________
__________ __________ __________ __________ __________
__________
__________
10
47 94 48 96 49 98
3 3 3 3 3 3
____________ ____________ ____________ ____________ ____________ ____________
__________
__________
__________
(e) For indices 1 and 2, estimate the order of convergence by taking the sequence of ratios of the error for N subintervals divided by the error for (2*N) subintervals and guessing the power of two that best approximates the limit of the sequence. (f) Estimate the order of accuracy for index 3 by taking the ve ratios of errors rk = ek /e2k for 49 k = 45, . . . , 49, take their geometric mean r = ( k=45 rk )1/5 , and guess the power of two that best approximates r. The reason that this case is dierent from the others is that the errors become near roundo and averaging is necessary to smooth out the resulting values. Geometric averaging is appropriate because the theoretical error curve is a straight line on a log-log plot and geometrical averaging is the same as arithmetic averaging of logs. The proofs you have seen about convergence of Gau quadrature rely on bounds on higher derivatives of the function. When bounds are not available, higher-order convergence might not be observed. Exercise 10: Consider the integral from Exercise 4
1
I=
0
log(x)dx =
/2.
(a) Use gl quad to ll in the following table. sqrt(-log(x)) Subingl_quad tervals index Error 10 20 40 80 1 1 1 1 ____________ ____________ ____________ ____________
Err ratio __________ __________ __________
What is the order of accuracy of the method using index=1? (b) Use gl quad to ll in the following table. sqrt(-log(x)) Subingl_quad tervals index Error 10 20 40 80 2 2 2 2 ____________ ____________ ____________ ____________
Err ratio __________ __________ __________
What is the order of accuracy of the method using index=2? (c) Use gl quad to ll in the following table.
11
sqrt(-log(x)) Subingl_quad tervals index Error 10 20 40 80 3 3 3 3 ____________ ____________ ____________ ____________
Err ratio __________ __________ __________
What is the order of accuracy of the method using index=3?
9
Adaptive quadrature
Our nal 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 atter. The advantage of this approach is that it minimizes the work necessary to compute a given integral. Ralston and Rabinowitz discuss adaptive quadrature starting on p. 126, and Atkinson discusses it starting on page 300. In this section, you will investigate a recursive algorithm for adaptively computing quadratures. 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. A good strategy for achieving a specied accuracy eciently is to attempt to uniformly distribute the error over each subinterval. If no interval is particularly bad, there is no obvious place to improve the estimate and if no interval is particularly good, no work has been wasted. Recall that, if an integration method has degree of precision p, then the local error on an integration interval of length h can be estimated by dividing the interval into two subintervals of length h/2 each. Denote by Qh the integral over the interval of length h and QL and QR the two estimates of the integral on the h/2 h/2 left and right subintervals. Then Qh = Q + Chp+2 f (p+1) () + h.o.t. QL + QR = Q + C(h/2)p+2 (f (p+1) ( L ) + f (p+1) ( R )) + h.o.t. h/2 h/2 (where C is a constant, , L , and R are appropriately chosen, and h.o.t. means higher order terms). Assuming that f (p+1) is roughly constant, f (p+1) () = f (p+1) ( L ) = f (p+1) ( R ) = f (p+1) , and that the higher order terms can be neglected yields Qh = Q + Chp+2 f (p+1) QL h/2 + QR h/2 = Q + C(h/2)
p+2
(4)
(p+1)
(2f
).
(5)
Solving Equation (4) for Q, plugging it into Equation (5), and dening the error as QL + QR Q yields h/2 h/2 the expression |QL + QR Qh | h/2 h/2 error estimate = (6) 2p+1 1 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. Exercise 11: 12
(a) Create a function m-le named adaptquad.m with the following code, and ll in the places marked ???. 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 % dont want xmid to appear twice in x x=[xleft xright(2:length(xright))]; Q= ??? ; errEst= ??? ; recursions=min(recursLeft,recursRight); end 13
Note: The input and output parameter recursions is not theoretically necessary, but is used to guard against innite recursion. The output vector x is not necessary, either, but will be used to show the eect of the adaptation. (b) Test adaptquad by applying it to the polynomial f5 (x) = 6x5 on the interval [0, 1], using tol=1.e-5 and recursions=50. Because Gauss-Legendre integration of index 3 is exact for this polynomial, the integral should equal 1 (i.e. error should be zero or roundo), and recursions=50, and the values of x should be three equally-spaced points in the interval. (c) Test adaptquad by applying it to the polynomial f6 (x) = 7x6 on the interval [0, 1], using tol=1.e-5 and recursions=50. Because Gauss-Legendre integration of index 3 is not exact for this polynomial, the integral should be close to 1. It turns out that a single set of renements is performed, so recursions=49, and the values of x should be ve equally-spaced points in the interval. The estimated and true errors should be extremely close. Please include the estimated error in your summary. (d) Test adaptquad by applying it to the Runge function on the interval [-5,5]. Use recursions=50. Recall that the exact answer is 2*atan(5). Fill in the following table adaptquad for Runge tol est. error exact error 1.e-3 __________ __________ 1.e-6 __________ __________ 1.e-9 __________ __________ You should nd that the estimated and exact errors are close in size, and smaller than tol. For the two more accurate cases, the estimated error is slightly larger than the exact error. As you can see, the estimated error is not so good for the case that tol=1.e-3. Exercise 12: Consider the following situation. A quadrature is being attempted with the call [Q,estErr,x,recursions]=adaptquad(func,0,1,tol,50); The estimated error is larger than tol at rst, so new calls with tol/2 are made for the intervals Ileft = [0, 0.5] and Iright = [0.5, 1]. Assume that the call for the interval Ileft satises the convergence criterion. Assume that the call for the interval Iright does not satisfy the convergence criterion, thus requiring two more calls to adaptquad. Assume that each of these calls satises the convergence criterion. What are the nal values of x, and recursions, after the adaptquad function has completed its work? The next few exercises will help you look a little more closely at the results of this recursive adaptive algorithm. Some of the points that will be made are listed below. You will see the advantage of adaptive algorithms. They save work over xed algorithms such as gl quad. You will see the pattern of the integration points take. They can be distributed in a very nonuniform fashion. You will see what happens when you try harder integrands. You will see some of the weaknesses of the algorithm.
14
In the following exercises, you will examine two functions that are more dicult to integrate. The rst is a scaled version of the Runge function, 1/(a + x2 ), where a = 105 The value of the integral on [-1,1] of the scaled Runge function is 1 1 2 1 dx = tan1 . 2 a a 1 a + x The scaled Runge function has a peak value of 1/a = 105 , so is much more strongly peaked than the unscaled Runge function. The second function is |x|, and the value of its integral over the interval [-1,1] is
1
|x|dx =
1
4 . 3
This function has a singularity in its derivatives at x = 0, thus invalidating the proof that the error estimator is reliable. A third function that is dicult to integrate is x0.99 . This function has an integrable singularity at x = 0 that is close to being nonintegrable. Exercise 13: (a) Write a function m-le named srunge.m for the scaled Runge function f (x) = 1/(a + x2 ) with a = 105 . (b) Evaluate the integral of srunge on the interval [-1,1] using the following call. [Q,estErr,x,recursions]=adaptquad(srunge,-1,1,1.e-10,50); What are the estimated and true errors? Is recursions larger than zero? (c) Use gl quad with index=3 on the scaled Runge function. Use trial-and-error to nd the value of numSubintervals that yields an error roughly comparable to estErr from adaptquad. How does this compare with length(x)-1, the number of subintervals that adaptquad used? (d) Plot the sizes of the subintervals that adaptquad used with the following command. 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. (e) What are the lengths of the largest and smallest intervals? Explain (one sentence) where you might expect to nd the smallest intervals for an arbitrary function. Exercise 14: (a) Approximate the integral of the function f (x) = |x| over the interval [-1,1] to a tolerance of 1.e-10. The exact value of this integral is 4/3. What are the estimated and true errors? (b) What is the returned value of recursions? It should be positive, indicating that the subinterval convergence criterion was always reached successfully. (c) Plot the subinterval sizes using the following command 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. 15
In the following exercise you will apply adaptquad to a function that is almost not integrable. You will see the benet of the recursions variable, whose reduction to 0 indicates that convergence was not achieved. Exercise 15: (a) Use adaptquad to approximate the integral
1 0
x0.99 dx = 100
to a tolerance of 1.e-10. Use recursions=50. Notice that the integration interval is [0,1]. What is the computed value of the integral? What are the estimated error and the true error? What is the value returned for recursions? (b) Do the same approximation starting with recursions=100. What are the computed value of the integral, the estimated error, the true error, and the returned value of recursions? (You may have to use the command set(0,RecursionLimit,200) in order that Matlab will allow sucient recursions.) Warning: This may take more than a minute on the computers in GSCC. (c) If the variable recursions were not used, this recursive function would never terminate because the convergence test would never be satised. (Actually, it would abort because there is a practical limit on the recursion depth.) The reason is that in the presence of the singularity, halving the interval does result in a reduction of error, but the reduction is half or less (c.f. Exercise 10), so it never passes the convergence test.
16