T he Filon − Simpson Code or : Numerical Integration of Highly Oscillatory Integrals R. Rosenfelder (PSI) June 1, 2006 Outline: 1. Introduction: Simpson, Newton, Gauss, Euler 2. Filon’s integration formula 3. New quadrature rules for a class of oscillatory integrals 4. Results 5. Summary
R. Rosenfelder (PSI) : T he Filon − Simpson Code 3 1. Introduction: Simpson, Newton, Gauss, Euler Numerical Integration: one-dimensional multidimensional functional ( aliter ...) ( totaliter aliter ...) ⇓ finite limits infinite limits singularities of integrand oscillatory integrands . . . fixed number of integration points adaptive (automatic) integration
R. Rosenfelder (PSI) : T he Filon − Simpson Code 4 Some integration rules: Trapezoidal rule ( h = x 1 − x 0 , f i ≡ f ( x i )) � x 1 − h 3 � � dx f ( x ) = h 12 f ′′ ( ξ ) , x 0 < ξ < x 1 f 0 + f 1 2 x 0 Extended trapezoidal rule ( h = x N − x 0 , x 0 < ξ < x N ) N � x N � 1 � − Nh 3 2 f 0 + f 1 + . . . + f N − 1 + 1 12 f ′′ ( ξ ) dx f ( x ) = h 2 f N x 0 Simpson’s rule � x 2 − h 5 � � dx f ( x ) = h 90 f (4) ( ξ ) f 0 + 4 f 1 + f 2 3 x 0
R. Rosenfelder (PSI) : T he Filon − Simpson Code 5 Extended Simpson’s rule � x 2 N � h dx f ( x ) = f 0 + 4 ( f 1 + f 3 + . . . + f 2 N − 1 ) 3 x 0 − Nh 5 � 90 f (4) ( ξ ) +2 ( f 2 + f 4 + . . . + f 2 N − 2 ) + f 2 N Newton-Cotes formulas (examples) � x 4 − 8 h 7 2 h � � 945 f (6) ( ξ ) closed type: dx f ( x ) = 7 f 0 + 32 f 1 + 12 f 2 + 32 f 3 + 7 f 4 45 x 0 � x 4 − 28 h 5 � � 4 h 90 f (4) ( ξ ) open type: dx f ( x ) = 2 f 1 − f 2 + 2 f 3 3 x 0
R. Rosenfelder (PSI) : T he Filon − Simpson Code 6 Gauss-Legendre integration � +1 N � dx f ( x ) = w i f ( x i ) + R N − 1 i =1 with 0 , w i = 2 [ P ′ N ( x i )] 2 P N ( x i ) = 1 − x 2 i 2 2 N +1 ( N !) 2 (2 N + 1)[(2 N )!] 3 f (2 N ) ( ξ ) , R N = − 1 < ξ < +1 Arbitrary interval: y = b + a + b − a 2 x 2 Different weight functions = ⇒ Gauss-Laguerre, Gauss-Hermite etc.
R. Rosenfelder (PSI) : T he Filon − Simpson Code 7 Euler & tanh sinh-integration from Borwein et al. : Experimentation in Mathematics: Computational Paths to Discovery (2004) p. 309 “Read Euler, read Euler ! He is the master of all of us” (Laplace) Euler-Maclaurin summation formula � b N f ( a + ih ) − h � dx f ( x ) = h 2 [ f ( a ) + f ( b )] a i =0 m h 2 j B 2 j � � � f (2 j − 1) ( b ) − f (2 j − 1) ( a ) − + R m (2 j )! j =1 − h 2 m +2 ( b − a ) B 2 m +2 f (2 m +2) ( ξ ) with R m = (2 m + 2)!
R. Rosenfelder (PSI) : T he Filon − Simpson Code 8 When f ( x ) and all its derivatives vanish at the endpoints a and b then the error of step-function approximation to integral is R m for all m → error goes to zero more rapidly than any power of h ! − Achieved by special transformation x = g ( t ) � +1 � + ∞ + ∞ � dt g ′ ( t ) f ( g ( t )) = h dx f ( x ) = w j f ( x j ) + R − 1 −∞ j = −∞ w j = g ′ ( jh ) with x j = g ( jh ) and For functions which are analytic in a strip | Im t | < d � e − πd/h � R = O Recommended transformation: x = tanh [ λ sinh t ] , t ∈ [ −∞ , + ∞ ] , λ > 0
R. Rosenfelder (PSI) : T he Filon − Simpson Code 9 Fast convergence also for | j | → ∞ since � − λ e | j | h � w j − → 2 λ exp Examples: ( λ = 1) � π/ 2 � 1 � 1 1 h dx sin x 0 dx ( − ln x ) 0 dx 2 √ x 0 1.0000 0.997 818 863 766 531 1.006 175 595 0 1.000 562 2 0.5000 0.999 999 556 735 862 0.999 999 700 1 0.999 999 5 0.2500 1.000 000 000 000 006 0.999 999 998 9 0.999 991 5 0.1250 0.999 999 999 999 999 0.999 999 990 6 0.999 978 2
R. Rosenfelder (PSI) : T he Filon − Simpson Code 10 2. Filon’s integration formula L. N. G. Filon: “On a quadrature formula for trigonometric integrals”, Proc. Royal Soc. Edinburgh 49 (1928), 38 – 47 � x 2 � � � dx f ( x ) cos xy = h α ( hy ) − f 0 sin( x 0 y ) + f 2 sin( x 2 y ) x 0 + β ( hy ) 1 � � f 0 cos( x 0 y ) + f 2 cos( x 2 y ) 2 � + γ ( hy ) f 1 cos( x 1 y ) + R with − 2 sin 2 Θ → 2Θ 3 Θ + sin 2Θ 1 Θ → 0 α (Θ) = − 45 + . . . 2Θ 2 Θ 3 � 1 + cos 2 Θ � − sin 2Θ → 2 Θ → 0 β (Θ) = 2 3 + . . . − Θ 2 Θ 3 � � − cos Θ + sin Θ → 4 Θ → 0 γ (Θ) = 4 − 3 + . . . Θ 2 Θ 3 ⇒ Simpson’s rule, analogous formula for � x 2 i.e. for y → 0 : = x 0 dx f ( x ) sin xy
R. Rosenfelder (PSI) : T he Filon − Simpson Code 11 How to derive this formula ? Write � x 2 dx f ( x ) cos xy ≃ w 0 f 0 + w 1 f 1 + w 2 f 2 x 0 and require that it is exact for polynomials x 0 , x 1 , x 2 . = ⇒ system of 3 equations for 3 unknowns � x 2 2 � dx x k cos( xy ) =: J k , w i x k = k = 0 , 1 , 2 , elementary integrals ! i x 0 i =0 Result: 1 w 0 = 2 h 2 [ x 1 x 2 J 0 − ( x 1 + x 2 ) J 1 + J 2 ] 1 w 1 = 2 h 2 [ − 2 x 0 x 2 J 0 + 4 x 1 J 1 − 2 J 2 ] 1 w 2 = 2 h 2 [ x 0 x 1 J 0 − ( x 0 + x 1 ) J 1 + J 2 ]
R. Rosenfelder (PSI) : T he Filon − Simpson Code 12 Asymptotic behaviour for y → ∞ : exact: � x 2 � x 2 � f ( x ) sin( xy ) x 2 − 1 dx f ′ ( x ) sin( xy ) � dx f ( x ) cos( xy ) = � y y x 0 x 0 x 0 � cos � � � 1 y →∞ − → f 2 sin( x 2 y ) − f 0 sin( x 0 y ) + O y 2 y Filon: � � � � y →∞ − → h α ( hy ) − f 0 sin( x 0 y ) + f 2 sin( x 2 y ) + . . . � �� � → 1 / ( hy ) is correct since the endpoints are included in quadrature formula !
R. Rosenfelder (PSI) : T he Filon − Simpson Code 13 3. New quadrature rules for a class of oscillatory integrals “In mathematics, I recognize true scientific value only in concrete truths, or to put it more pointedly, only in mathematical (computations)” (Kronecker) Motivation: Worldline variational approach to QFT: describe relativistic particles by their worldlines x µ ( t ) , t = proper time and not by fields Φ( x ) S t ∼ � dt � dt ′ ˙ x ( t ) A ( t − t ′ ) ˙ x ( t ′ ) Quadratic trial action + Feynman-Jensen variational principle � exp( − ∆ S ) � ≥ exp( −� ∆ S � ) = ⇒ pair of non-linear variational equations for � ∞ � Eσ � 1 + const . 1 δV δµ 2 ( σ ) sin 2 A ( E ) = dσ “profile function” E 2 2 0 � ∞ sin 2 ( Eσ/ 2) 1 µ 2 ( σ ) ∼ dE “pseudotime” E 2 A ( E ) 0
R. Rosenfelder (PSI) : T he Filon − Simpson Code 14 V [ µ 2 ] ∼ � S I � is the interaction term averaged over the trial action, specific for the field theory under consideration, δV [ µ 2 ] → const . σ → 0 − δµ 2 ( σ ) σ 2+ r the ”force” ( r = 0 for super-renormalizable, r = 1 for renormalizable theories) �� � 2 � µ 2 ( σ = t − t ′ ) ∼ x ( t ) − x ( t ′ ) mean square displacement � ∞ dµ 2 ( σ ) 1 sin( Eσ ) ∼ dE mean ”velocity” (also needed for QED) dσ A ( E ) E 0 Large- E, σ limit of A ( E ) and µ 2 ( σ ) may be obtained analytically = ⇒ restriction to finite integration limits , add the analytically calculated asymptotic contribution to the numerical result
R. Rosenfelder (PSI) : T he Filon − Simpson Code 15 Therefore consider � b dx f ( x ) sin( xy ) I 1 [ f ]( a, b, y ) : = xy a � b dx f ( x ) 4 sin 2 ( xy/ 2) I 2 [ f ]( a, b, y ) : = x 2 y 2 a i.e. � b I j [ f ]( a, b, y ) = dx f ( x ) O j ( xy ) , O j (0) = 1 a Numerical evaluation of Fourier integrals with infinite upper limits is much more demanding: see NAG routine D01ASF based on strategy in: Piessens et al.: QUADPACK, A Subroutine Package for Automatic Integration , Springer (1983) requires a delicate extrapolation procedure, not suitable for a numerical solu- tion of the variational equations
R. Rosenfelder (PSI) : T he Filon − Simpson Code 16 “As long as right methods are used, quadrature of highly-oscillatory integrals is very accurate and affordable !” (A. Iserles) strategy: Choose N points x ( j ) in the interval [ a, b ], two of them identical i with the endpoints and require that the integral over x k O j ( xy ) is exact. For simplification: N = 2, equally spaced points x i = a + ih = ⇒ as in Filon’s case � � 1 w ( j ) x 1 x 2 J ( j ) − ( x 1 + x 2 ) J ( j ) + J ( j ) = 0 0 1 2 2 h 2 1 � � w ( j ) − 2 x 0 x 2 J ( j ) + 4 x 1 J ( j ) − 2 J ( j ) = 1 0 1 2 2 h 2 1 � � w ( j ) x 0 x 1 J ( j ) − ( x 0 + x 1 ) J ( j ) + J ( j ) = 2 0 1 2 2 h 2 but now moments involve special functions � � 1 J ( j ) F ( j ) ( by ) − F ( j ) = ( ay ) k k k y k +1 � z dt t k O j ( t ) with F ( j ) ( z ) = k 0
R. Rosenfelder (PSI) : T he Filon − Simpson Code 17 We have F (1) Si(z) , F (1) 1 (z) = 1 − cos z , F (1) ( z ) = 2 (z) = sin z − z cos z 0 � � Si(z) − 1 − cos z F (2) , F (2) ( z ) = 2 ( z ) = 2 [ γ + ln z − Ci(z) ] , 0 1 z F (2) ( z ) = 2 [ z − sin z ] 2 with γ = 0 . 5772156640 . . . (Euler’s constant) and � z dt sin t Si(z) : = Sine integral t 0 � z dt cos t − 1 Ci(z) : = γ + ln z + Cosine integral t 0 y → 0: obtain Simpson weights = ⇒ new quadrature rules may be called Filon-Simpson rules Divide [ a, b ] in N intervals and apply FS-rules in each interval = ⇒ Extended Filon-Simpson quadrature rules
Recommend
More recommend
