Generalized Techniques in Numerical Integration Richard M. - PowerPoint PPT Presentation

Generalized Techniques in Numerical Integration Richard M. Slevinsky and Hassan Safouhi Mathematical Section Campus Saint-Jean, University of Alberta hassan.safouhi@ualberta.ca http://www.csj.ualberta.ca/safouhi Approximation and Extrapolation of Convergent and Divergent Sequences and Series CIRM Luminy September 28 – October 2, 2009 Generalized Techniques in Numerical Integration. – p. 1/29

The Plan Introduction A Mathematical Tool Formulae for Higher Order Derivatives PART I – Ongoing Research The general Idea Example of Applications and Numerical Results PART II – Almost Completed An Algorithm for the G (1) Transformation n Computing the Incomplete Bessel Functions Generalized Techniques in Numerical Integration. – p. 2/29

Introduction The Euler series arising from integrating the Euler integral by parts: � ∞ � ∞ � ∞ e − t e − x e − t ( − 1) l l ! x l + ( − 1) n n ! t d t = t n +1 d t x x x l =0 � ∞ e − x ( − 1) l l ! ∼ x → ∞ . x l , x l =0 Integration by parts by x d x led to: � d � λ � sin( x ) � � ∞ � ∞ g ( x ) ( − 1) λ x λ g ( x ) j λ ( x ) d x = d x x d x x 0 0 � � d � λ � sin( x ) � � ∞ � x λ − 1 g ( x ) = ( − 1) λ x d x. x d x x 0 Generalized Techniques in Numerical Integration. – p. 3/29

Introduction �� � � d � λ − 1 � sin( x ) � ∞ ( − 1) λ � � x λ − 1 g ( x ) = � x d x x � 0 � d � � � � d � λ − 1 � sin( x ) � � ∞ x λ − 1 g ( x ) + ( − 1) λ − 1 x d x x d x x d x x 0 �� � d � � d � l � � λ − 1 − l � sin( x ) � ∞ λ − 1 � � x λ − 1 g ( x ) ( − 1) λ + l = � � x d x x d x x l =0 0 � d � λ � � � sin( x ) � � ∞ x λ − 1 g ( x ) + x d x. x d x x 0 Leading at: �� d �� � λ � � ∞ � ∞ 1 x λ − 1 g ( x ) g ( x ) j λ ( v x ) d x = sin( v x ) d x. v λ +1 x d x 0 0 Generalized Techniques in Numerical Integration. – p. 4/29

Introduction Semi-infinite spherical Bessel integrals in molecular integrals: ˆ n k n + 1 2 [ R γ ( s, x )] � 2 ( z ) = z n ( n + j )! 1 ˆ g ( x ) = x n x , k n + 1 [ γ ( s, x )] n γ e z (2 z ) j j !( n − j )! j =0 � (1 − s ) ζ 2 i + sζ 2 j + s (1 − s ) x 2 . s ∈ [0 , 1] . γ ( s, x ) = 40 0.04 0.04 Integrand with spherical Bessel Integrand with spherical Bessel Integrand with sine function 0.03 0.03 30 0.02 0.02 20 0.01 0.01 0 0 10 -0.01 -0.01 0 -0.02 -0.02 -0.03 -0.03 -10 -0.04 -0.04 -20 �� 0 0 2 2 4 4 6 6 8 8 10 10 0 0.5 1 1.5 2 2.5 3 3.5 4 � b f ( x ) d x ? How can we use this technique for any a Generalized Techniques in Numerical Integration. – p. 5/29

Higher Order Derivatives Let us determine the k th derivatives of G 1 ( x ) = x 3 f ( x 2 ) : � d � G 1 ( x ) = 3 x 2 f ( x 2 ) + 2 x 4 f ′ ( x 2 ) . d x � d � 2 G 1 ( x ) = 6 x f ( x 2 ) + (6 x 3 + 8 x 3 ) f ′ ( x 2 ) + 4 x 5 f ′′ ( x 2 ) . d x How about this: � d � d � � k ( x − 3 G 1 ( x )) = 2 k f ( k ) ( x 2 ) . ( x − 3 G 1 ( x )) = 2 f ′ ( x 2 ) = ⇒ x d x x d x � � k d For G 2 ( x ) = x 2 f (ln( x )) : ( x − 2 G 2 ( x )) = f ( k ) (ln( x )) . x − 1 d x � d � k � � i d ( x − n G ( x )) for G ( x ) in terms of Can we express x m d x d x i = 0 , 1 , ..., k ? Generalized Techniques in Numerical Integration. – p. 6/29

Higher Order Derivatives The Slevinsky-Safouhi formulae [Slevinsky and Safouhi, 2009]: � � k d ( x − n G ( x )) Theorem Let G ( x ) be k th differentiable with x m d x well-defined. The Slevinsky-Safouhi formula I for ( α, β, m, n ) is given by: � � k � � i k � d d ( x − β G ( x ))= A i k x n − β + i ( m +1) − k ( α +1) ( x − n G ( x )) , x α d x x m d x i =0 with coefficients [ N = ( n − β − ( k − 1)( α + 1)) ]:   1 i = k for  A i N A 0 k = i = 0 , k > 0 for k − 1   k − 1 + A i − 1 ( N + i ( m + 1)) A i 0 < i < k. for k − 1 i ( − 1) i − j ( n − β + j ( m + 1) − ( k − 1)( α + 1)) k,α +1 � A i k = , m � = − 1 . ( m + 1) i j ! ( i − j )! j =0 The Slevinsky-Safouhi formula II: ( α, β, m, n ) = (0 , 0 , 1 , 0) . Generalized Techniques in Numerical Integration. – p. 7/29

Part I The Generalized S n Transformation & The Staircase Algorithm Generalized Techniques in Numerical Integration. – p. 8/29

The generalized S n � b Let f ( x ) be integrable on [ a, b ] , i.e. a f ( x ) d x exists. We write: � b � b f ( x ) d x = G 0 ( x ) H 0 ( x ) w ( x ) d x, a a for some weight function w ( x ) , whose choice depends on f ( x ) . � b If f ( x ) ∈ C n [ a, b ] , then a f ( x ) d x has the equivalent representation, which we obtain after n integration by parts by w ( x ) d x : � � b � b n − 1 � b � � f ( x ) d x = G l ( x ) H l + 1 ( x ) + G n ( x ) H n ( x ) w ( x ) d x � a a a l = 0 � � l � � − l where: d d G l ( x ) = ( − 1) l g ( x ) H l ( x ) = h ( x ) . and w ( x ) d x w ( x ) d x Generalized Techniques in Numerical Integration. – p. 9/29

The Staircase Algorithm � b f ( x ) d x take the following form: Approximations to a � x 0 For a < x 0 < b , initialize: S 0 = G 0 ( x ) H 0 ( x ) w ( x ) d x, a � � b � b b � � G 0 ( x ) H 0 ( x ) w ( x ) d x = G 0 ( x ) H 1 ( x ) + G 1 ( x ) H 1 ( x ) w ( x ) d x. � x 0 x 0 x 0 � � x 1 b � � S 1 = S 0 + G 0 ( x ) H 1 ( x ) + G 1 ( x ) H 1 ( x ) w ( x ) d x, x 0 < x 1 < b. � x 0 x 0 For the sequence { x l } n l =1 satisfying a < x l − 1 < x l < b , define: � � x l b � � S l = S l − 1 + G l − 1 ( x ) H l ( x ) + G l ( x ) H l ( x ) w ( x ) d x. � x l − 1 x l − 1 � b f ( x ) d x form the sequence { S l } n The approximations to l =0 . a Generalized Techniques in Numerical Integration. – p. 10/29

Bessel Integral The integral that follows appeared in Numerical Recipes: � b x I 1 = x 2 + 1 J 0 ( x ) d x = K 0 (1) . 0 1 By choosing w ( x ) = x , we have G 0 ( x ) = x 2 + 1 and H 0 ( x ) = J 0 ( x ) . 2 l l ! H l ( x ) = x l J l ( x ) . → ֒ G l ( x ) = and ( x 2 + 1) l +1 The integral then has the equivalent representations: � � ∞ n − 1 − 2 l l ! x l + 1 � x n + 1 � ∞ + 2 n n ! � I 1 = ( x 2 + 1 ) l + 1 J l + 1 ( x ) ( x 2 + 1 ) n + 1 J n ( x ) d x � a 0 l = 0 All the boundary terms vanish and consequently: � ∞ � ∞ x n +1 x x 2 + 1 J 0 ( x ) d x = 2 n n ! ( x 2 + 1) n +1 J n ( x ) d x = K 0 (1) . 0 0 Generalized Techniques in Numerical Integration. – p. 11/29

Bessel Integral - Results Table 1: I 1 = 0 . 421024438240708 . x l = 2 π ( l + 1) . ¯ ¯ l S l l S l 0 0.414 193 276 771 795 7 0.421 024 438 053 642 1 0.421 696 746 593 657 8 0.421 024 438 183 915 2 0.421 072 353 906 909 9 0.421 024 438 241 771 3 0.421 020 653 966 770 10 0.421 024 438 241 364 4 0.421 023 974 243 269 11 0.421 024 438 240 707 5 0.421 024 464 857 404 12 0.421 024 438 240 700 6 0.421 024 443 256 394 13 0.421 024 438 240 708 Generalized Techniques in Numerical Integration. – p. 12/29

Fresnel Integrals The integrals are given by: � ∞ � ∞ ˜ sin( v x 2 ) d x cos( v x 2 ) d x. I 2 ( a, v ) = I 2 ( a, v ) = and a a By choosing w ( x ) = x , we have G 0 ( x ) = 1 x and H 0 ( x ) = sin( v x 2 ) . H l ( x ) = sin( v x 2 − lπ/ 2) (2 l )! → ֒ G l ( x ) = . and 2 l l ! x 2 l +1 (2 v ) l The integral I 2 ( a, v ) then has the equivalent representations: � � ∞ sin( vx 2 − ( l + 1 ) π � sin( vx 2 − n π n − 1 � ) 2 ) − 2 ( 2l )! ( 2n )! � 2 + d x � ( 4v ) l + 1 l ! x 2l + 1 ( 4v ) n n ! x 2n � a l = 0 a Generalized Techniques in Numerical Integration. – p. 13/29

Fresnel Integrals - Results � Table 2: I 2 (0 , 1) = . 626657068657750 . x l = 2 π ( l + 1) /v . ¯ ¯ l S l l S l 0 0.629 878 864 869 732 7 0.626 657 068 644 466 1 0.627 294 419 199 049 8 0.626 657 068 657 872 2 0.626 651 451 723 302 9 0.626 657 068 657 790 3 0.626 655 488 072 699 10 0.626 657 068 657 749 4 0.626 657 083 038 776 11 0.626 657 068 657 749 5 0.626 657 073 115 469 12 0.626 657 068 657 750 6 0.626 657 068 616 796 The integral ˜ I 2 ( a, v ) then has the equivalent representations: � � ∞ cos( v x 2 − ( l + 1 ) π � cos( vx 2 − n π n − 1 � ) 2 ) − 2 ( 2l )! ( 2n )! � 2 + d x � ( 4v ) l + 1 l ! x 2l + 1 ( 4v ) n n ! x 2n � a l = 0 a Generalized Techniques in Numerical Integration. – p. 14/29

The Twisted Tail The integral is proposed in the book " The SIAM 100-Digit Challenge ": � 1 � ∞ � � t − 1 cos t − 1 ln( t ) cos( x e x ) d x, I 3 = ( x = − ln( t )) . d t = 0 0 1 w ( x ) = (1 + x ) e x ֒ (1 + x ) e x and H 0 ( x ) = cos( x e x ) . → G 0 ( x ) = � � l � � − d 1 x e x − lπ ֒ → G l ( x )= (1 + x ) e x and H l ( x ) = cos . (1 + x ) e x d x 2  The general form of G l ( x ) is:  p 0 ( x ) = 1      p 1 ( x ) = 2 + x  e − ( l +1) x 9 + 8 x + 2 x 2 G l ( x ) = (1 + x ) 2 l +1 p l ( x ) p 2 ( x ) =   64 + 79 x + 36 x 2 + 6 x 3  p 3 ( x ) =    . . .  . . . . . . Generalized Techniques in Numerical Integration. – p. 15/29