A comparative study of extrapolation methods, sequence - PowerPoint PPT Presentation

A comparative study of extrapolation methods, sequence transformations and steepest descent methods Computing infinite-range integrals Richard M. Slevinsky and Hassan Safouhi Mathematical Section Campus Saint-Jean, University of Alberta hsafouhi@ualberta.ca SC2011 – International Conference on Scientific Computing October 10 – 14, 2011

Oscillatory integrals: A numerical challenge − 3 x 10 8 0.15 6 0.1 4 0.05 2 0 0 − 0.05 − 2 − 0.1 − 4 − 0.15 − 6 − 0.2 − 8 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Classical quadrature deteriorates rapidly as the oscillations become strong. New methods / New oscillatory quadrature methods were developed.

The plan The three most popular methods : Extrapolation methods Sequence transformations Numerical steepest descent Applications & Comparison : We computed four integrals using the three methods We introduced some refinements to the algorithms We performed comparisons with regard to efficiency Conclusion

Extrapolation methods Let f ( x ) satisfy: m ∞ α i � � p k ( x ) f ( k ) ( x ) p k ( x ) ∼ x k f ( x ) = where as x → ∞ . x i k =1 i =0 Then [ Levin & Sidi 1981 ]: � ∞ m − 1 ∞ β k , i � � x k +1 f ( k ) ( x ) f ( t ) d t ∼ x → ∞ . as x i x i =0 k =0 � ∞ Let D ( m ) represent approximations to f ( t ) d t [ Levin & Sidi ]: n 0 � x l m − 1 n − 1 ¯ β k , i D ( m ) � x k +1 f ( k ) ( x l ) � = f ( t ) d t + , n l x i 0 l k =0 i =0 where { x l } mn +1 is an increasing sequence. l =0 In general, m is equal or can be reduced to 1 or 2.

Extrapolation methods D (2) When m = 2, we use the ¯ approximation [ Sidi 1982 ]: n � x n − 1 ¯ β 1 , i D (2) ¯ � = F ( x l ) + x 2 l f ′ ( x l ) F ( x ) = f ( t ) d t , where n x i 0 l i =0 where { x l } n +1 l =0 are the successive positive zeros of f ( x ). To compute D (1) D (2) or ¯ n , we use the W algorithm [ Sidi 1982 ]: n = M (1) = M (2) WD (1) D (2) n W ¯ n . or n n N (1) N (2) n n M ( j ) and N ( j ) for n , j = 1 , 2 , . . . , are computed recursively by: n n F ( x j ) 1 M ( j ) N ( j ) = and = 0 0 x 2 x 2 j f ′ ( x j ) j f ′ ( x j ) M ( j +1) − M ( j ) N ( j +1) − N ( j ) M ( j ) N ( j ) n − 1 n − 1 n − 1 n − 1 = and = . n n x − 1 j + n − x − 1 x − 1 j + n − x − 1 j j

Sequence transformations Let us consider I ( λ ): � ∞ ∞ � I ( λ ) = f ( x ; λ ) d x ∼ a k ( λ ) . 0 k =0 S [ a ] − S n [ a ] ∼ R n [ a ] ⇒ S [ a ] ≈ S n [ a ] + R n [ a ] . Consider the case where the remainder is given by: ∞ c j � S [ a ] − S n [ a ] ∼ ω n n → ∞ . as ( n + β ) j j =0 A Levin transformation [ Levin 1973 ]: � ( n + β + j ) k − 1 k � k S n + j [ a ] � ( − 1) j ( n + β + k ) k − 1 ω n + j j j =0 L ( n ) k ( β ) = . � ( n + β + j ) k − 1 k � k 1 � ( − 1) j ( n + β + k ) k − 1 j ω n + j j =0

Sequence transformations A recursive algorithm [ Fessler et al. 1983 ] to compute L ( n ) k ( β ): 1 For n = 0 , 1 , . . . , set: = S n [ a ] = 1 P ( n ) Q ( n ) and . 0 0 ω n ω n 2 For n = 0 , 1 , . . . , k = 1 , 2 , . . . , compute P ( n ) and Q ( n ) from: k k � k − 2 β + n � β + n + k − 1 U ( n ) = U ( n +1) U ( n ) − k − 1 , k k − 1 β + n + k β + n + k where the U ( n ) stand for either P ( n ) or Q ( n ) k . k k 3 For all n and k , set: = P ( n ) L ( n ) k . k Q ( n ) k We choose ω n = a n , which gives rise to the t ( n ) k ( β ) transformation.

Numerical steepest descent Huybrechs and Vandewalle 2006 . Consider the integral: � b e i w g ( x ) = e − w ℑ g ( x ) e i w ℜ g ( x ) . f ( x ) e i w g ( x ) d x , I = a The steepest descent is based on: 1 e i w g ( x ) decays exponentially fast if ℑ g ( x ) > 0. 2 e i w g ( x ) does not oscillate if ℜ g ( x ) is fixed. 3 The value of I does not depend on the exact path taken (Cauchy theorem). A new path is defined at a [ Huybrechs & Vandewalle 2006 ]: g ( h a ( κ )) = g ( a ) + κ i κ ≥ 0 . with The integral is then equivalent to: � ∞ f ( a + κ i ) e − w κ i d κ. I = F ( a ) − F ( b ) where F ( t ) = e i w g ( a ) 0

Numerical steepest descent 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 − 0.2 − 0.2 − 0.4 − 0.4 − 0.6 − 0.6 � − 0.8 − 0.8 0 5 10 15 20 25 30 0 5 10 15 20 25 30

The first integral � ∞ e − x x + β d x = − e β Ei ( − β ) . I 1 ( β ) = 0 No substitution is required to apply Gauss-Laguerre quadrature. The integrand satisfies a first order linear differential equation: x + β x + β + 1 f ′ f 1 ( x ) = − 1 ( x ) . The integral has the asymptotic expansion: ∞ I 1 ( β ) ∼ 1 k ! � as β → ∞ . β ( − β ) k k =0

The first integral Table: Numerical valuation of I 1 ( β ). Error WD (1) Error LT (0) Error GL n β n n n n n 0.03 124 .87(-03) 17 .22(-12) 16 .32(-01) 0.10 124 .25(-05) 17 .84(-12) 17 .11(-02) 0.30 124 .16(-09) 17 .71(-14) 17 .14(-05) 1.00 124 .13(-13) 12 .56(-15) 21 .18(-07) 3.00 34 .36(-14) 8 .13(-14) 20 .21(-12) 4.00 34 .40(-14) 7 .22(-14) 19 .40(-15) 5.00 22 .65(-15) 6 .31(-14) 18 .16(-15) 10.00 15 .61(-15) 3 .30(-14) 16 .30(-15) 30.00 9 .17(-14) 3 .47(-14) 13 .00( 00) 100.00 6 .18(-15) 2 .21(-14) 9 .18(-15) Calculation time WD (1) Calculation time LT (0) n n = 0 . 49 and Calculation time GL n = 0 . 075. Calculation time GL n

Refinement – The first integral For the sequence transformation : As the values of the governing parameter(s) tend to 0 + , we resort to using a series representation of the integrals: � ∞ � ( − β ) k � I 1 ( β ) = − e β β → 0 + . C + ln β + as k · k ! k =1 For the Steepest descent : � ∞ � x n � ∞ f ( x ) e − x d x = e − x d x + e − x n f ( x + x n ) e − x d x 0 0 0 � x n � x i n e − x d x = f ( x ) e − x d x . � Compute 0 x i − 1 i =1 n is determined by: � x n � x n − 1 � � � � � < 1 � � � � f ( x ) e − x d x f ( x ) e − x d x � . � � � � 2 � � � � x n − 1 x n − 2 � �

Refinement – The first integral Table: Numerical valuation of I 1 ( β ) with the refinement. Error WD (1) Error LT (0) Error GL n β n n n n n 0.03 70 .10(-12) 17 .22(-12) 7 .15(-15) 0.10 53 .18(-14) 17 .84(-12) 9 .00( 00) 0.30 41 .13(-14) 17 .71(-14) 12 .18(-15) 1.00 22 .13(-14) 12 .56(-15) 17 .00( 00) 3.00 12 .17(-14) 8 .13(-14) 28 .47(-14) 4.00 12 .17(-14) 7 .22(-14) 19 .40(-15) 5.00 12 .28(-14) 6 .31(-14) 18 .16(-15) 10.00 8 .33(-14) 3 .30(-14) 16 .30(-15) 30.00 6 .47(-14) 3 .47(-14) 13 .00( 00) 100.00 5 .19(-14) 2 .21(-14) 9 .18(-15) Calculation time WD (1) Calculation time LT (0) n n = 2 . 2 and Calculation time GL n = 0 . 27. Calculation time GL n

The second integral � ∞ √ sin( b x ) d x = π � a � a � I 2 ( a , b ) = sin bJ 1 (2 a b ) . 2 x 0 Asymptotic expansion as a b → ∞ : √ � π a 1 � I 2 ( a , b ) ∼ cos(2 a b − 3 π/ 4) √ 2 b � 2 a b ∞ ( − 1) k Γ(3 / 2 + 2 k ) � × (16 a b ) k (2 k )! Γ(3 / 2 − 2 k ) k =0 √ ∞ � ( − 1) k − sin(2 a b − 3 π/ 4) Γ(5 / 2 + 2 k ) � √ . (16 a b ) k (2 k + 1)! Γ(1 / 2 − 2 k ) 4 a b k =0

The second integral � a Splitting the integration interval with respect to x 0 = b : � x 0 � ∞ � a � a � � I 2 ( a , b ) = sin sin( b x ) d x + sin sin( b x ) d x x x 0 x 0 � sin( a x ) � ∞ � ∞ � b � a � = sin d x + sin sin( b x ) d x x 2 x x x − 1 x 0 0 � e i a x �� ∞ �� ∞ � � b � a � � e i b x d x = Im sin + Im sin . x 2 d x x x x − 1 x 0 0 The substitutions x new = i a ( x − 1 − x old ) and x new = i b ( x 0 − x old ) 0 lead to: � ∞ � � � e i a x − 1 � 0 e − x i b I 2 ( a , b ) = sin + i x / a ) 2 d x Im x − 1 ( x − 1 a + i x / a 0 0 0 � i � ∞ � a � � e i b x 0 e − x d x + sin . Im b x 0 + i x / b 0

The second integral Table: Numerical valuation of I 2 ( a , b ). Error W ¯ D (2) Error LT (0) Error GL n a b n n n n n 1.0 1.0 124 .19(-09) 16 .86(-15) 23 .29(-08) 2.0 1.0 124 .60(-11) 17 .12(-14) 24 .34(-10) 2.0 2.0 124 .72(-12) 17 .17(-14) 24 .33(-12) 3.0 1.0 124 .12(-11) 18 .27(-15) 23 .10(-11) 3.0 2.0 124 .75(-14) 17 .55(-15) 25 .17(-14) 3.0 3.0 117 .68(-14) 17 .51(-15) 18 .13(-15) 10.0 1.0 124 .29(-14) 19 .16(-14) 20 .44(-15) 100.0 1.0 124 .76(-14) 10 .91(-01) 8 .21(-14) 10.0 10.0 124 .59(-14) 7 .79(-02) 8 .20(-14) 100.0 10.0 69 .20(-14) 8 .26( 01) 5 .16(-14) Calculation time W ¯ D (2) Calculation time LT (0) n n = 0 . 096 and Calculation time GL n = 0 . 011. Calculation time GL n