Numerical Integration Sanzheng Qiao Department of Computing and - PowerPoint PPT Presentation

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Numerical Integration Sanzheng Qiao Department of Computing and Software McMaster University August, 2012

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Outline Introduction 1 Rectangle Rule 2 Trapezoid Rule 3 Error Estimates 4 Simpson’s Rule 5 Richardson’s Extrapolation 6 Adaptive Quadrature 7 2D Quadrature 8 Software Packages 9

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Introduction Another (better) term: quadrature Problem: Given finite number of function values f ( x i ) , x i ∈ [ a , b ] or the function f ( x ) can be evaluated at any x ∈ [ a , b ] , calculate � b I ( f ) = f ( x ) dx . a

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Introduction (cont.) Partition a = x 1 < x 2 < · · · < x n + 1 = b , and denote h i = x i + 1 − x i . Then � x i + 1 n I ( f ) = I i I i = f ( x ) dx � x i i = 1 Quadrature rule: Approximation of I i Composite quadrature rule: Approximation of I ( f ) as a sum of I i

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Rectangle rule We use piecewise constant (degree zero polynomial) to approximate f ( x ) . In each interval [ x i , x i + 1 ] , f ( x ) is evaluated at the midpoint y i = x i + x i + 1 i = 1 , ..., n , , 2 then the rectangle (quadrature) rule is: I i ≈ h i f ( y i ) x x i i+1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Rectangle rule (cont.) The composite rectangle rule is: n R ( f ) = h i f ( y i ) � i = 1 A weighted sum of function values. Often the major computation is the evaluation of the function. Thus the complexity is measured by the number of function evaluations. In the rectangle rule Function evaluations: n . Evaluated at midpoints y i .

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Trapezoid rule We use piecewise linear interpolation (degree one polynomial) to approximate f ( x ) . In each interval [ x i , x i + 1 ] the function is evaluated at the endpoints: f ( x i ) + f ( x i + 1 ) I i ≈ h i 2 x x i i+1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Trapezoid rule (cont.) Composite trapezoid rule: n f ( x i ) + f ( x i + 1 ) T ( f ) h i � = 2 i = 1 h 1 2 f ( x 1 ) + h 1 + h 2 f ( x 2 ) + ... = 2 + h n − 1 + h n f ( x n ) + h n + 1 f ( x n + 1 ) 2 2 A weighted sum of function values. In the trapezoid rule, Function evaluations: n + 1. Evaluated at endpoints x i

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Error in rectangle rule Taylor expansion f ( x ) about the midpoint y i = ( x i + x i + 1 ) / 2: ( x − y i ) p ∞ f ( p ) ( y i ) . f ( x ) = f ( y i ) + � p ! p = 1 Integrate the both sides and note that � x i + 1 h p + 1 � i ( x − y i ) p dx = even p ( p + 1 ) 2 p x i 0 odd p

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Error in rectangle rule (cont.) Then � x i + 1 f ( x ) dx x i h i f ( y i ) + 1 1 i f iv ( y i ) + · · · 24 h 3 i f ′′ ( y i ) + 1920 h 5 = When h i is small, the error I ( f ) − R ( f ) n n 1 1 i f iv ( y i ) h 3 i f ′′ ( y i ) + h 5 � � ≈ 24 1920 i = 1 i = 1 For equal spacing, h i = h , we have n n I ( f ) − R ( f ) ≈ h 3 h 5 f iv ( y i ) f ′′ ( y i ) + � � 24 1920 i = 1 i = 1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Error in trapezoid rule In order to make the error in the trapezoid rule comparable with that in the rectangle rule, we expand f ( x ) at the midpoint y i . Substituting x = x i and x = x i + 1 in the Taylor expansion, we have ( − 1 ) p h p ∞ 2 p p ! f ( p ) ( y i ) i f ( x i ) f ( y i ) + � = p = 1 h p ∞ 2 p p ! f ( p ) ( y i ) i f ( x i + 1 ) f ( y i ) + � = p = 1 Thus f ( x i ) + f ( x i + 1 ) = f ( y i ) + 1 1 i f iv ( y i ) + · · · 8 h 2 i f ′′ ( y i ) + 384 h 4 2

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Error in trapezoid rule (cont.) Recall that in the case of rectangle rule, we had � x i + 1 f ( x ) dx x i h i f ( y i ) + 1 1 i f iv ( y i ) + · · · 24 h 3 i f ′′ ( y i ) + 1920 h 5 = Combining the above two equations, we have � x i + 1 f ( x ) dx x i f ( x i ) + f ( x i + 1 ) − 1 1 i f iv ( y i ) + · · · h i 12 h 3 i f ′′ ( y i ) − 480 h 5 = 2 Then the error is n n I ( f ) − T ( f ) ≈ − 1 1 i f iv ( y i ) + · · · � h 3 i f ′′ ( y i ) − � h 5 12 480 i = 1 i = 1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Remarks Compare n n I ( f ) − R ( f ) ≈ h 3 h 5 f iv ( y i ) f ′′ ( y i ) + � � 24 1920 i = 1 i = 1 and n n I ( f ) − T ( f ) ≈ − 1 1 i f iv ( y i ) + · · · � h 3 i f ′′ ( y i ) − � h 5 12 480 i = 1 i = 1 Usually rectangle rule (degree zero approximation) is more accurate than trapezoid rule (degree one approximation). Surprised?

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Remarks Observe n n I ( f ) − R ( f ) ≈ h 3 h 5 f iv ( y i ) f ′′ ( y i ) + � � 24 1920 i = 1 i = 1 and n n I ( f ) − T ( f ) ≈ − 1 1 i f iv ( y i ) + · · · h 3 i f ′′ ( y i ) − h 5 � � 12 480 i = 1 i = 1 Using I ( f ) − R ( f ) ≈ 1 3 ( T ( f ) − R ( f )) , we can estimate the error in R ( f ) using T ( f ) and R ( f ) . Similarly, I ( f ) − T ( f ) ≈ 2 3 ( R ( f ) − T ( f )) can be used to estimate the error in T ( f ) . (But they are approximations, it is possible that R ( f ) − T ( f ) = 0 whereas I ( f ) − R ( f ) � = 0).

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Remarks Observe n n I ( f ) − R ( f ) ≈ h 3 h 5 f iv ( y i ) f ′′ ( y i ) + � � 24 1920 i = 1 i = 1 When each h i is cut in half, I ( f ) − R 1 2 ( f ) ≈ 1 4 ( I ( f ) − R ( f )) . (Why?) Similarly for the trapezoid rule. Doubling the number of panels in either the rectangle rule or the trapezoid rule, it can be expected to roughly quadruple the accuracy. This can be used to estimate the error as well as improving the accuracy. (How?)

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Example Compute π � 2 sin ( x ) dx 0 using the trapezoid rule. m QCTrap(sin, 0.000, 1.571, m) error 3 0.9480594489685199 5.2e-2 5 0.9871158009727753 1.3e-2 9 0.9967851718861696 3.2e-3 17 0.9991966804850722 8.0e-4 where m is the number of points, that is, m − 1 is the number of intervals.

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Simpson’s rule Recall the rectangle rule n n R ( f ) = I ( f ) − 1 1 i f iv ( y i ) + · · · h 3 i f ′′ ( y i ) − h 5 � � 24 1920 i = 1 i = 1 and the trapezoid rule n n T ( f ) = I ( f ) + 1 1 i f iv ( y i ) + · · · h 3 i f ′′ ( y i ) + h 5 � � 12 480 i = 1 i = 1 Combining the above two equations (canceling the O ( h 3 i ) term), we get a more accurate method (Simpson’s rule): 2 3 R ( f ) + 1 S ( f ) 3 T ( f ) = n 1 i f iv ( y i ) + · · · I ( f ) + h 5 � = 2880 i = 1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Simpson’s rule (cont.) Simpson’s rule: 3 h i f ( x i + x i + 1 f ( x i ) + f ( x i + 1 ) I i = 2 ) + 1 3 h i 2 2 Composite Simpson’s rule: n f ( x i ) + 4 f ( x i + x i + 1 1 � � S ( f ) = 6 h i ) + f ( x i + 1 ) � 2 i = 1 Function evaluations: 2 n + 1 Error n 1 i f iv ( y i ) + · · · I ( f ) − S ( f ) = − h 5 � 2880 i = 1

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary Remarks Simpson’s rule can also be derived by using piecewise quadratic (degree two) approximation. Actually, Simpson’s rule is exact for cubic function (one extra order of accuracy), since the error term involves the fourth derivatives. Doubling the number of panels in Simpson’s rule can be expected to reduce the error by roughly the factor of 1 / 16.

Intro Rectangle Trapezoid Error Simpson’s Extrapolation Adaptive 2D Software Summary A general technique: Richardson’s extrapolation Idea: Combining two approximations (e.g., R ( f ) and T ( f ) ) which have similar error terms to achieve a more accurate approximation (e.g., S ( f ) ). Example. Combining S ( f ) and S 1 2 ( f ) to obtain an approximation which has error of order h 7 i . This gives the Romberg quadrature. Question What are the weights? 15 S 1 2 ( f ) − 1 15 S ( f ) Answer: 16