Numerical integration (a.k.a quadrature formulas or quadrature rules) Quadrature rules are used to approximate integrals of functions that we are not able to compute exactly. Given g : [ c , d ] → R , the most common quadrature rules look like � d k +1 � g ( x ) dx ≃ ω i g ( x i ) c i =1 where: x 1 , x 2 , · · · , x k +1 are the quadrature “ points” or “nodes” of the rule and ω 1 , ω 2 , · · · , ω k +1 are called quadrature “weights” Definition: The order of precision of a quadrature rule is the maximum degree of the polynomials which are integrated exactly by the rule. Among the numerous quadrature rules, we shall see the so-called interpolatory rules. December 20, 2019 1 / 15
Interpolatory quadrature rules The function g is approximated by its Lagrange interpolant Π k ( x ) (of degree ≤ k =) with respect to the given nodes x 1 , x 2 , · · · , x k +1 , and then the integral of the polynomial is computed exactly. Then � d � d � d k +1 � g ( x ) dx ≃ Π k ( x ) dx = g ( x i ) L i ( x ) dx c c c i =1 � d k +1 � = g ( x i ) L i ( x ) dx c i =1 � �� � ω i The order of precision of an interpolator formula will be at least k (if g ∈ P k , it is integrated exactly) December 20, 2019 2 / 15
Error analysis for interpolatory quadrature rules Bounds for the quadrature error are derived by the bounds for the interpolation error (see the slides on Lagrange interpolation): x ∈ [ c , d ] | g ( x ) − Π k ( x ) | ≤ ( d − c ) k +1 x ∈ [ c , d ] | g ( k +1) ( x ) | max max ( ∗ ) ( k + 1)! Thus, � d � d � d � d � � � � � � � � g ( x ) dx − Π k ( x ) dx � = ( g ( x ) − Π k ( x )) dx � ≤ | g ( x ) − Π k ( x ) | dx � � � � c c c c � d ≤ max [ c , d ] | g ( x ) − Π k ( x ) | dx = ( d − c ) max [ c , d ] | g ( x ) − Π k ( x ) | (1) c ≤ ( d − c ) k +2 x ∈ [ c , d ] | g ( k +1) ( x ) | max ( k + 1)! (Observe that if g ∈ P k , then g ( k +1) ( x ) ≡ 0. Hence, the quadrature error is = 0) This estimate is obtained using the generic bound for the interpolation error. Sharper estimates can be obtained by analysing each rule directly, as we shall see. December 20, 2019 3 / 15
Use of quadrature rules It should be clear by now that if we want a good approximation of an integral we have to use properly the rules in order to make the error as smaller as we want. Exactly like we did for Lagrange interpolation, we will construct piecewise integration rules, also called composite integration rules. Given f : [ a , b ] → R (smooth enough), subdivide [ a , b ] in N subintervals, for simplicity of notation all equal. We have then a uniform subdivision of [ a , b ] into intervals of length h = ( b − a ) / N : I 1 = [ x 1 , x 2 ] , · · · , I j = [ x j , x j +1 ] , · · · , I N = [ x N , x N +1 ]. In each subinterval we approximate f with a Lagrange interpolant polynomial of degree k . Thus, � b � x j +1 � x j +1 N N � � f ( x ) dx = f ( x ) dx ≃ Π k ( x ) dx a x j x j j =1 j =1 Let us see some examples. December 20, 2019 4 / 15
Composite midpoint rule x M =midpoint of the interval I j : x M = ( x j + x j +1 ) / 2 j j f ( x ) | [ a , b ] ≃ f 0 ( x ) piecewise constant function given by f 0 ( x ) | I j = f ( x M j ) j = 1 , 2 , · · · , N � b � b � x j +1 N N � � f ( x M f ( x M f ( x ) dx ≃ f 0 ( x ) dx = j ) dx = h j ) a a x j j =1 j =1 December 20, 2019 5 / 15
Composite midpoint rule: pseudocode on the blackboard December 20, 2019 6 / 15
Composite trapezoidal rule f ( x ) | [ a , b ] ≃ f 1 ( x ) piecewise linear function which, on each interval I j , is the Lagrange interpolant of degree ≤ 1 with respect to the endpoints of I j � b � b � x j +1 N N Π 1 ( x ) dx = h � � f ( x ) dx ≃ f 1 ( x ) dx = [ f ( x j ) + f ( x j +1 )] 2 a a x j j =1 j =1 December 20, 2019 7 / 15
Exercice on composite midpoint and trapezoidal rules f Exercise " and the function Cx ) the given g. : × = following E- I ] mesh 1 on : , I h = - Iz . +13 ×¢=y " I ×2= Xs = - . ' far ) of f. compute the approximation by a , midpoint rule composite • rule trapezoidal composite • December 20, 2019 8 / 15
Exercice on composite midpoint and trapezoidal rules f Exercise A ) " and the function the given g : x = following E- I ] mesh 1 on : , xin xin Xiu f- h - - 13 1/4=1 Xie I Xz t Fg Xs - = - ' far ) of f. compute the approximation by a , midpoint rule composite [ • it yourself rule do trapezoidal composite ← • ( quadrature points ) XF 43,27-0,4723 the midpoints 3 are = and rule the is : " ) h ( if f. fkn tf ) ) G) dx= Cx 'd ) tfcxjn Zz ( Ft 'T ) 6¥ - = 91 ¥ fCx7 ) f CAY and " f ) C x 's o since - = = = = December 20, 2019 9 / 15
Composite Simpson rule f ( x ) | [ a , b ] ≃ f 2 ( x ) piecewise quadratic function which, on each interval I j , is the Lagrange interpolant of degree ≤ 2 with respect to the endpoints and the midpoint of I j � b � b N f 2 ( x ) dx = h � [ f ( x j ) + 4 f ( x M f ( x ) dx ≃ j ) + f ( x j +1 )] 6 a a j =1 December 20, 2019 10 / 15
Composite Midpoint rule: error bound � x j +1 N � � � � ( f ( x ) − f ( x M � ERR = j )) dx � � x j j =1 � �� � E j The first possibility is to use (1), which gives in our case: � � � x j +1 � � � ≤ h 2 ′ ( x ) | = O ( h 2 ) � � E j = f ( x ) − Pi 0 ( x ) dx x ∈ [ x j , x j +1 ] | f max � � � x j Then the global error ERR is the sum of the error in each subinterval: N � [ a , b ] | f ′ ( x ) | h 2 N = max [ a , b ] | f ′ ( x ) | ( b − a ) h = O ( h ) ERR ≤ | E j | ≤ max j =1 (we used Nh = b − a and max [ I j ] | f ′ ( x ) | ≤ max [ a , b ] | f ′ ( x ) | ). December 20, 2019 11 / 15
Composite Midpoint rule: error bound (sharp) � x j +1 N � � � � ( f ( x ) − f ( x M � ERR = j )) dx � � x j j =1 � �� � E j Second possibility: use in each I j the Taylor expansion centered in x M j ( x − x M j ) 2 f ( x ) = f ( x M j ) + ( x − x M j ) f ′ ( x M f ′′ ( z ) ( for z between x and x M j ) + j ) 2! � x j +1 � x j +1 j ) 2 ( x − x M ( x − x M j ) f ′ ( x M f ′′ ( z ) dx E j = j ) dx + 2! x j x j � �� � = 0 � x j +1 max [ I j ] | f ′′ ( x ) | ≤ 1 j ) 2 dx = [ I j ] | f ′′ ( x ) | h 3 ( x − x M 2 max 24 x j December 20, 2019 12 / 15
Composite Midpoint rule: error bound (sharp) Again, the global error is the sum of the error in each subinterval. N | E j | ≤ max [ a , b ] | f ′′ ( x ) | h 3 N = max [ a , b ] | f ′′ ( x ) | � ( b − a ) h 2 ERR ≤ 24 24 j =1 (we used Nh = b − a and max [ I j ] | f ′′ ( x ) | ≤ max [ a , b ] | f ′′ ( x ) | ). We see that the error is zero if f ′′ ≡ 0 in each I j , i.e., if f is a piecewise polynomial of degree 1. Hence, the order of precision of the midpoint rule is actually 1 (even though we are projecting onto piecewise constants!) To summarize: with C = max [ a , b ] | f ′′ ( x ) | ERR ≤ C h 2 ( b − a ) 24 We see that ERR → 0 for h → 0 quadratically with h (halving h reduces the error by 4). December 20, 2019 13 / 15
Hints on Gaussian rules * NOT FOR THE EXAM * Conclusions The order of precision of an interpolatory quadrature rule using n nodes is at least n − 1 (that is, the rule integrates exactly polynomials of degree up to n − 1). December 20, 2019 14 / 15
Hints on Gaussian rules * NOT FOR THE EXAM * Conclusions The order of precision of an interpolatory quadrature rule using n nodes is at least n − 1 (that is, the rule integrates exactly polynomials of degree up to n − 1). If the n nodes are Gauss points (special points that are defined as roots of Legendre polynomials) the order of precision is higher: precisely, 2 n − 1. December 20, 2019 14 / 15
Hints on Gaussian rules * NOT FOR THE EXAM * Conclusions The order of precision of an interpolatory quadrature rule using n nodes is at least n − 1 (that is, the rule integrates exactly polynomials of degree up to n − 1). If the n nodes are Gauss points (special points that are defined as roots of Legendre polynomials) the order of precision is higher: precisely, 2 n − 1. Operatively: Gauss points are computed in the open interval ] − 1 , 1[: � 1 n = 1 g ( x ) dx ≈ 2 g (0) − 1 � 1 � � � � − 1 + 1 n = 2 g ( x ) dx ≈ g √ + g √ 3 3 − 1 √ � √ � � � � 1 g ( x ) dx ≈ 5 3 + 8 9 g (0) + 5 3 √ √ n = 3 9 g − 9 g 5 5 − 1 n = 4 ... see the reference books December 20, 2019 14 / 15
Hints on Gaussian rules * NOT FOR THE EXAM * To compute Gauss points on a generic interval [ a , b ] and use them for � b evaluating g ( x ) dx is simple. Consider the map a F : [ − 1 , 1] → [ a , b ] , that is, ˆ x ∈ [ − 1 , 1] → x = F (ˆ x ) ∈ [ a , b ] December 20, 2019 15 / 15
Hints on Gaussian rules * NOT FOR THE EXAM * To compute Gauss points on a generic interval [ a , b ] and use them for � b evaluating g ( x ) dx is simple. Consider the map a F : [ − 1 , 1] → [ a , b ] , that is, ˆ x ∈ [ − 1 , 1] → x = F (ˆ x ) ∈ [ a , b ] It is easy to check that the map is linear, given by x = b − a x + a + b 2 ˆ 2 December 20, 2019 15 / 15
Recommend
More recommend
