Chap 9 Numerical Integration

Chap 9 Numerical Integration - FIRST YEAR CALCULUS W W L...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: FIRST YEAR CALCULUS W W L CHEN c W W L Chen, 1987, 2005. This chapter originates from material used by the author at Imperial College, University of London, between 1981 and 1990. It is available free to all individuals, on the understanding that it is not to be used for financial gain, and may be downloaded and/or photocopied, with or without permission from the author. However, this document may not be kept on any information storage and retrieval system without permission from the author, unless such system is not accessible to any individuals other than its owners. Chapter 9 NUMERICAL INTEGRATION 9.1. Introduction In Chapters 7 and 8, we have discussed some analytic techniques for evaluating integrals. However, many integrals that arise in science and engineering resist attack by even the most sophisticated analytic techniques. In such instances, we may have to accept a rather poor and perhaps even not entirely satisfactory second best, and attempt to make reasonable approximations by numerical techniques. 9.2. The Trapezium Rule Suppose that we wish to evaluate an integral B f (x) dx, A where the function f (x) is finite and continuous in the closed interval [A, B ]. If we draw the curve y = f (x), then the value of the integral is the same as the area bounded by the curve y = f (x) and the lines y = 0, x = A and x = B (the reader should draw a diagram). A first, and rather crude, approximation to the integral is to take the area of the trapezium with vertices at the points (A, 0), (B, 0), (A, f (A)) and (B, f (B )). In other words, we take the approximation B (1) A f (x) dx ≈ 1 (B − A)(f (A) + f (B )). 2 In practice, however, we take more points than just A and B . Consider the dissection A = x0 < x1 < . . . < xn = B Chapter 9 : Numerical Integration page 1 of 10 First Year Calculus c W W L Chen, 1987, 2005 of the interval [A, B ]. Clearly we have B n xi f (x) dx = A i=1 xi−1 f (x) dx. Suppose now that we make a similar approximation as (1) in each of the subintervals, so that for every i = 1, . . . , n, we have xi xi−1 f (x) dx ≈ 1 (xi − xi−1 )(f (xi−1 ) + f (xi )). 2 Then we have the approximation B (2) A f (x) dx ≈ 1 2 n (xi − xi−1 )(f (xi−1 ) + f (xi )). i=1 Suppose further that the lengths of all the subintervals are the same, so that xi − xi−1 = h = Then (2) becomes B B−A n for every i = 1, . . . , n. f (x) dx ≈ A h (f (x0 ) + 2f (x1 ) + 2f (x2 ) + . . . + 2f (xn−1 ) + f (xn )). 2 This is called the Trapezium rule for n intervals. Example 9.2.1. We wish to estimate the value of log 2 by a Trapezium rule approximation to the integral 2 1 1 dx. x Then f (x) = 1/x in the interval [1, 2]. If we take h = 1/2, then we have x f (x) and so 2 1 1 1 3 2 2 3 2 1 2 1 1 dx ≈ x 4 1+ 41 + 32 = 0.7083 (4dp). If we take h = 1/4, then we have x f (x) and so 2 1 Chapter 9 : Numerical Integration 1 1 5 4 4 5 3 2 2 3 7 4 4 7 2 1 2 1 1 dx ≈ x 8 1+ 8481 +++ 5372 = 0.6970 (4dp). page 2 of 10 First Year Calculus c W W L Chen, 1987, 2005 9.3. The Midpoint Rule This method is fundamentally similar to the Trapezium rule. Suppose that we wish to evaluate an integral B f (x) dx, A where the function f (x) is finite and continuous in the closed interval [A, B ]. Consider the point C = 1 (A + B ), the midpoint in the interval [A, B ] (the reader should draw a 2 diagram). We take the approximation B (3) A f (x) dx ≈ (B − A)f A+B 2 . Suppose that we divide the interval [A, B ] into n subintervals by the dissection A = x0 < x1 < . . . < xn = B, and make a similar approximation as (3) in each of the subintervals, so that for every i = 1, . . . , n, we have xi xi−1 f (x) dx ≈ (xi − xi−1 )f xi−1 + xi 2 . Then we have the approximation B n (4) A f (x) dx ≈ i=1 (xi − xi−1 )f xi−1 + xi 2 . Suppose further that the lengths of all the subintervals are the same, so that xi − xi−1 = h = Then (4) becomes B n B−A n for every i = 1, . . . , n. f (x) dx ≈ h A i=1 f xi−1 + xi 2 . This is called the Midpoint rule for n intervals. Example 9.3.1. We wish to estimate the value of log 2 by a Midpoint rule approximation to the integral 2 1 1 dx. x Then f (x) = 1/x in the interval [1, 2]. If we take h = 1/2, then we have x f (x) and so 2 1 Chapter 9 : Numerical Integration (1) 5 4 4 5 (3) 2 7 4 4 7 (2) 1 1 dx ≈ x 2 44 + 57 = 0.6857 (4dp). page 3 of 10 First Year Calculus c W W L Chen, 1987, 2005 If we take h = 1/4, then we have x f (x) and so 2 1 (1) 9 8 8 9 (5) 4 11 8 8 11 (3) 2 13 8 8 13 (7) 4 15 8 8 15 (2) 1 1 dx ≈ x 4 8 8 8 8 + + + 9 11 13 15 = 0.6912 (4dp). 9.4. Simpson’s Rule Suppose that we wish to evaluate an integral B f (x) dx, A where the function f (x) is finite and continuous in the closed interval [A, B ]. In both the Trapezium rule and the Midpoint rule, a crude approximation to the area under the curve is obtained by replacing the curve between x = A and x = B by a straight line segment; in other words, a polynomial of degree 1. A natural extension of this idea is to replace the curve by a parabola; in other words, a polynomial of degree 2, passing through the points (A, f (A)), (B, f (B )) and (C, f (C )), where C = 1 (A + B ). 2 Consider first the simple case A = −H and B = H , so that C = 0. We wish to fit a parabola p(x) = αx2 + βx + γ through these points. Then αH 2 − βH + γ = f (−H ), γ = f (0), 2 αH + βH + γ = f (H ), so that α= f (−H ) − 2f (0) + f (H ) , 2H 2 f (H ) − f (−H ) β= , 2H γ = f (0). We now take the approximation H −H H f (x) dx ≈ (αx2 + βx + γ ) dx = −H 2 H αH 3 + 2γH = (f (−H ) + 4f (0) + f (H )). 3 3 In general, if we wish to use this approximation over the interval [A, B ], the above becomes B f (x) dx ≈ A Chapter 9 : Numerical Integration B−A 6 f (A) + 4f A+B 2 + f (B ) . page 4 of 10 First Year Calculus c W W L Chen, 1987, 2005 This is called Simpson’s rule with 3 ordinates. Similar to the Trapezium rule and the Midpoint rule, we may apply Simpson’s rule on subintervals of the interval [A, B ]. Suppose that we divide the interval [A, B ] into n subintervals by the dissection A = x0 < x1 < . . . < xn = B, where n is even and xi − xi−1 = h = B−A n for every i = 1, . . . , n. Applying Simpson’s rule to the interval [x2j −2 , x2j ], we have x2j x2j −2 f (x) dx ≈ x2j − x2j −2 (f (x2j −2 ) + 4f (x2j −1 ) + f (x2j )), 6 n−1 n−1 so that B f (x) dx ≈ A h f (x0 ) + f (xn ) + 4 3 f (xi ) + 2 i=1 i odd i=1 i even f (xi ) . This is called Simpson’s rule for (n + 1) ordinates, where n is even. Note that the coefficients for f (x0 ), f (x1 ), f (x2 ), . . . , f (xn ) are respectively 1, 4, 2, 4, 2, 4, 2, . . . , 4, 2, 4, 1. Example 9.4.1. We wish to estimate the value of log 2 by a Simpson rule approximation to the integral 2 1 1 dx. x Then f (x) = 1/x in the interval [1, 2]. If we take h = 1/2, hence 3 ordinates, then we have x f (x) and so 2 1 1 1 3 2 2 3 2 1 2 1 1 dx ≈ x 6 1+ 81 + 32 = 0.6944 (4dp). If we take h = 1/4, hence 5 ordinates, then we have x f (x) and so 2 1 1 1 5 4 4 5 3 2 2 3 7 4 4 7 2 1 2 1 1 dx ≈ x 12 1+ 16 4 16 1 ++ + 5 3 7 2 = 0.6933 (4dp). page 5 of 10 Chapter 9 : Numerical Integration First Year Calculus c W W L Chen, 1987, 2005 9.5. Truncation Errors In this section, we state without proof some results concerning the errors that inevitably occur when we apply the Trapezium rule, Midpoint rule or Simpson’s rule. As far as numerical integration is concerned, such error analysis is more important than the estimates that are given by the rules. The study of these questions forms part of numerical analysis. For the Trapezium rule and Midpoint rule, we have the following two results. PROPOSITION 9A. Suppose that the function f (x) is finite and continuous in the closed interval [A, B ], and that the second derivative f (x) exists for every x ∈ (A, B ). Suppose further that the Trapezium rule, applied to the dissection A = x0 < x1 < . . . < xn = B of [A, B ] into n subintervals, where xi − xi1 = h = gives rise to the error B B−A n for every i = 1, . . . , n, Tn = A f (x) dx − h (f (x0 ) + 2f (x1 ) + 2f (x2 ) + . . . + 2f (xn−1 ) + f (xn )). 2 Then |Tn | ≤ where K = max |f (x)|. x∈[A,B ] K (B − A)3 K (B − A)h2 , = 2 12n 12 Furthermore, if f (x) does not change sign in the interval [A, B ], then Tn has the opposite sign to the sign of f (x) in this interval. PROPOSITION 9B. Suppose that the function f (x) is finite and continuous in the closed interval [A, B ], and that the second derivative f (x) exists for every x ∈ (A, B ). Suppose further that the Midpoint rule, applied to the dissection A = x0 < x1 < . . . < xn = B of [A, B ] into n subintervals, where xi − xi−1 = h = gives rise to the error B n B−A n for every i = 1, . . . , n, Mn = A f (x) dx − h i=1 f xi−1 + xi 2 . Then |Mn | ≤ Chapter 9 : Numerical Integration K (B − A)3 K (B − A)h2 , = 2 24n 24 page 6 of 10 First Year Calculus c W W L Chen, 1987, 2005 where K = max |f (x)|. x∈[A,B ] Furthermore, if f (x) does not change sign in the interval [A, B ], then Mn has the same sign as the sign of f (x) in this interval. Example 9.5.1. In our Trapezium and Midpoint rule approximation to log 2, we have used the function f (x) = 1/x in the interval [1, 2]. Note that f (x) = 2/x3 > 0 in this interval. It follows that Tn < 0 and Mn > 0. This means that our Trapezium rule estimates are over-estimates, and our Midpoint rule estimates are under-estimates. The corresponding result for Simpson’s rule is somewhat different. PROPOSITION 9C. Suppose that the function f (x) is finite and continuous in the closed interval [A, B ], and that the fourth derivative f (x) exists for every x ∈ (A, B ). Suppose further that Simpson’s rule, applied to the dissection A = x0 < x1 < . . . < xn = B of [A, B ] into n subintervals, where n is even and xi − xi−1 = h = gives rise to the error B B−A n for every i = 1, . . . , n, n−1 n−1 Sn = A f (x) dx − h f (x0 ) + f (xn ) + 4 3 f (xi ) + 2 i=1 i odd i=1 i even f (xi ) . Then |Sn | ≤ where L = max |f x∈[A,B ] L(B − A)5 L(B − A)h4 , = 4 180n 180 (x)|. Furthermore, if f (x) does not change sign in the interval [A, B ], then Sn has the opposite sign to the sign of f (x) in this interval. Example 9.5.2. In our Simpson rule approximation to log 2, we have used the function f (x) = 1/x in the interval [1, 2]. Note that f (x) = 24/x5 > 0 in this interval. It follows that Sn < 0. This means that our estimates are over-estimates. 9.6. Richardson Extrapolation Throughout, we assume that the function f (x) is continuous in the closed interval [A, B ], and write B I= A Chapter 9 : Numerical Integration f (x) dx. page 7 of 10 First Year Calculus c W W L Chen, 1987, 2005 Consider first of all the Trapezium rule. Suppose that f (x) does not change sign in the interval [A, B ]. Suppose further that T (h) denotes the Trapezium rule approximation to I with a given h. Then in view of Proposition 9A, we have I − T (h) ≈ CT (B − A)h2 , where CT is a constant, so that I − T (h) ≈ CT (B − A). h2 Repeating the same argument on the Trapezium rule approximation to I with h/2, we have I − T (h/2) ≈ CT (B − A). (h/2)2 It follows that I − T (h) I − T (h/2) ≈ , h2 (h/2)2 so that I − T (h) ≈ 4(I − T (h/2)), whence I≈ 4T (h/2) − T (h) . 3 Example 9.6.1. Recall our estimates of log 2 in Example 9.2.1. We have T (1/2) = 0.7083 (4dp) Hence I≈ 4(0.6970) − (0.7083) ≈ 0.6932. 3 and T (1/4) = 0.6970 (4dp). Consider next the Midpoint rule. Suppose that f (x) does not change sign in the interval [A, B ]. Suppose further that M (h) denotes the Midpoint rule approximation to I with a given h. Then in view of Proposition 9B, we have I − M (h) ≈ CM (B − A)h2 , where CM is a constant, so that I − M (h) ≈ CM (B − A). h2 Repeating the same argument on the Midpoint rule approximation to I with h/2, we have I − M (h/2) ≈ CM (B − A). (h/2)2 It follows that I − M (h) I − M (h/2) ≈ , h2 (h/2)2 Chapter 9 : Numerical Integration page 8 of 10 First Year Calculus c W W L Chen, 1987, 2005 so that I − M (h) ≈ 4(I − M (h/2)), whence I≈ 4M (h/2) − M (h) . 3 Example 9.6.2. Recall our estimates of log 2 in Example 9.3.1. We have M (1/2) = 0.6857 (4dp) Hence I≈ 4(0.6912) − (0.6857) ≈ 0.6930. 3 and M (1/4) = 0.6912 (4dp). Consider finally Simpson’s rule. Suppose that f (x) does not change sign in the interval [A, B ]. Suppose further that S (h) denotes the Simpson rule approximation to I with a given h. Then in view of Proposition 9C, we have I − S (h) ≈ CS (B − A)h4 , where CS is a constant, so that I − S (h) ≈ CS (B − A). h4 Repeating the same argument on the Simpson rule approximation to I with h/2, we have I − S (h/2) ≈ CS (B − A). (h/2)4 It follows that I − S (h) I − S (h/2) ≈ , h4 (h/2)4 so that I − S (h) ≈ 16(I − S (h/2)), whence I≈ 16S (h/2) − S (h) . 15 Example 9.6.3. Recall our estimates of log 2 in Example 9.4.1. We have S (1/2) = 0.6944 (4dp) Hence I≈ 16(0.6933) − (0.6944) ≈ 0.6932. 15 and S (1/4) = 0.6933 (4dp). Remark. It is worth noting that log 2 = 0.6931 (4dp). Chapter 9 : Numerical Integration page 9 of 10 First Year Calculus c W W L Chen, 1987, 2005 Problems for Chapter 9 1. Consider the integral a) b) c) d) e) f) g) h) i) 1 dx. 2 0 1+x Find the Trapezium rule approximation with 2 intervals. Find the Trapezium rule approximation with 4 intervals. Find the Midpoint rule approximation with 2 intervals. Find the Midpoint rule approximation with 4 intervals. Discuss whether the estimates in (a)–(d) are over-estimates or under-estimates. Justify your assertions. Find the Simpson rule approximation with 2 intervals. Find the Simpson rule approximation with 4 intervals. Use Richardson extrapolation on your results in (a)–(d). What number are we approximating? 1 Chapter 9 : Numerical Integration page 10 of 10 ...
View Full Document

This note was uploaded on 02/01/2009 for the course MATH 3423341 taught by Professor Staff during the Spring '09 term at UCSD.

Ask a homework question - tutors are online