Quadrature of highly oscillatory integrals: the role of (complex) - PowerPoint PPT Presentation

Quadrature of highly oscillatory integrals: the role of (complex) orthogonal polynomials Alfredo Dea˜ no Optimal Point Configurations and Orthogonal Polynomials 2017 CIEM, Castro Urdiales. April 18-22, 2017 Joint work with Daan Huybrechs (KU Leuven, Belgium) and Arieh Iserles (University of Cambridge, UK) Alfredo Dea˜ no Highosc and OPs

Outline Classical Gaussian quadrature and connection to OPs Highly oscillatory integrals Asymptotic analysis vs quadrature Numerical steepest descent Complex Gaussian quadrature Alfredo Dea˜ no Highosc and OPs

Gaussian quadrature Standard quadrature Q [ f ] of a 1 D integral has the form � b n � I [ f ] = f ( x ) w ( x )d x ≈ b k f ( x k ) = Q [ f ] , a k =1 where w ( x ) is a positive weight function , with finite moments x k are the quadrature nodes b k are the quadrature weights General question: How to choose the nodes and weights? General answer: It depends on what you want to do! Alfredo Dea˜ no Highosc and OPs

Gaussian quadrature I Some standard approaches: Set equispaced nodes h = b − a x k = a + kh, k = 0 , 1 , . . . , n, , n construct the Lagrange interpolating polynomial and integrate. This gives Newton-Cotes rules (trapezoidal, Simpson...) Use ideas from approximation theory. For example, maximise the degree of exactness: I [ p ] = Q [ p ] , p ∈ P m , for m as large as possible. Here P m = { p ( x ) polynomial , deg( P ) ≤ m } . Alfredo Dea˜ no Highosc and OPs

Gaussian quadrature II Theorem Given an integer k with 0 ≤ k ≤ n , the quadrature rule has degree of exactness equal to d = n − 1 + k if and only if The quadrature is interpolatory 1 The nodal polynomial 2 n � p n ( x ) = ( x − x k ) k =1 satisfies � b p n ( x ) q ( x ) w ( x )d x = 0 , q ∈ P k − 1 a This is maximised when k = n , then we have Gaussian quadrature. Alfredo Dea˜ no Highosc and OPs

Gaussian quadrature and OPs In Gaussian quadrature, one has to work with OPs p n ( x ) : � b p n ( x ) x k w ( x )d x = 0 , k = 0 , 1 , 2 , . . . , n − 1 , a where w ( x ) is positive and has finite moments. Classical cases: nodes and weights well understood and tabulated. Non-classical cases: how to compute nodes and weights? Linear Algebra methods (Golub–Welsch) Modified moments Discretization methods Asymptotic methods W. Gautschi. Orthogonal polynomials. Computation and Approximation. OUP, 2004. Alfredo Dea˜ no Highosc and OPs

Highly oscillatory integrals One problem: approximate the value of oscillatory integrals like � b f ( x )e i ωx d x, I [ f ] = a where ω ≫ 1 . Other examples are � sin( ωx ) � b � b � f ( x ) H (1 , 2) I [ f ] = f ( x ) d x, I [ f ] = ( ωx )d x, ν cos( ωx ) a a or more in general � b f ( x )e i ωg ( x ) d x, I [ f ] = a where f ( x ) and g ( x ) are assumed to be smooth. Alfredo Dea˜ no Highosc and OPs

Highly oscillatory integrals We can try Gaussian quadrature: consider � 1 � 1 1 1 2 + x e i ωx d x, 2 + x e i ωx 2 d x. I 1 [ f ] = I 2 [ f ] = − 1 − 1 Figure: Errors in log 10 scale in approximating I 1 [ f ] (left) and I 2 [ f ] with Gaussian quadrature for 0 ≤ ω ≤ 100 and n = 4 , 8 , 12 , . . . , 32 . Alfredo Dea˜ no Highosc and OPs

Gaussian quadrature vs asymptotics If we apply integration by parts, we obtain � 1 f ( x )e i ωx d x I [ f ] = − 1 s 1 � f ( k ) (1)e i ω − f ( k ) ( − 1)e − i ω � � + O ( ω − s − 1 ) = − ( − i ω ) k +1 k =0 as ω → ∞ . Some observations: This is an asymptotic expansion, hence divergent in general. However, it is bound to be very good for large ω ! It only uses information at the endpoints x = ± 1 . Alfredo Dea˜ no Highosc and OPs

Oscillatory integrals 1 0.5 0 −0.5 −1 0 1 2 3 4 5 6 7 8 9 10 1 0.5 0 −0.5 −1 0 1 2 3 4 5 6 7 8 9 10 Figure: Plot of cos(10 − x ) sin ωx , ω = 10 (top) and ω = 1000 (bottom). Alfredo Dea˜ no Highosc and OPs

Stationary points 1 0.5 0 −0.5 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1 0.5 0 −0.5 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Figure: Plot of cos(100 x 2 ) (top) and cos(100 x 3 ) (bottom). Alfredo Dea˜ no Highosc and OPs

Other methods I Filon-type methods: interpolate only the non–oscillatory part of the integrand. Interpolate where? At the points given by large ω asymptotics! Construct p ( x ) of degree 2 s + 1 that satisfies p ( i ) ( − 1) = f ( i ) ( − 1) , p ( i ) (1) = f ( i ) (1) , i = 0 , 1 , . . . , s, and then � 1 p ( x )e i ωx d x. Q [ f ] = − 1 We try the same example as before: � 1 � 1 1 1 2 + x e i ωx 2 d x. 2 + x e i ωx d x, I 1 [ f ] = I 2 [ f ] = − 1 − 1 Alfredo Dea˜ no Highosc and OPs

Other methods II Figure: Errors in log 10 scale in approximating I 1 [ f ] (left) and I 2 [ f ] with the Filon method 0 ≤ ω ≤ 100 and s = 0 , 1 , . . . , 6 . Alfredo Dea˜ no Highosc and OPs

Other methods III Extended Filon methods (recent work of J. Gao, A. Iserles): add extra internal nodes to improve results for small/moderate ω . Filon–Jacobi, using zeros of Jacobi polynomial P ( s +1 ,s +1) ( x ) . n Filon–Clenshaw–Curtis, with zeros of Chebyshev polynomials ( 1 2 , 1 2 ) ( x ) (Dom´ ınguez, Graham, Smyshlyaev 2011). P n Adaptive Filon. Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent A different alternative is the following: take � 1 f ( x )e i ωg ( x ) d x I [ f ] = − 1 and deform [ − 1 , 1] following paths of steepest descent of g ( x ) in C . For example, for g ( x ) = x : Im z C − 1 + i P i P 1 + i P Re z − 1 1 � 1 − 1 f ( x )e i ωx d x . Figure: Steepest descent paths for the oscillatory integral Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent We choose P > 0 and we write � 1 f ( x )e i ωx d x I [ f ] = − 1 � 1 � P f ( − 1 + i p )e i ω ( − 1+i p ) d p + f ( x + i P )e i ω ( x +i P ) d x = i 0 − 1 � P f (1 + i p )e i ω (1+i p ) d p − i 0 Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent We choose P > 0 and we write � 1 f ( x )e i ωx d x I [ f ] = − 1 � P = ie − i ω f ( − 1 + i p )e − ωp d p 0 � 1 � P + e − ωP f ( x + i P )e i ωx d x − ie i ω f (1 + i p )e − ωp d p. − 1 0 High oscillation turns into exponential decay! Then, we can discretise the path integrals using (for example) Gauss–Laguerre quadrature. Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent More complicated g ( x ) functions can be used, including stationary points (Asheim, Huybrechs 2010, 2013). Advantages: optimal asymptotic order in ω . Drawbacks: expensive computation and analyticity requirements. Alfredo Dea˜ no Highosc and OPs

Complex Gaussian quadrature An intriguing alternative: consider � 1 f ( x )e i ωx d x, I [ f ] = − 1 and apply Gaussian quadrature directly , with weight w ( x ) = e i ωx . Large ω asymptotic analysis in D., Huybrechs, Iserles 2015. Also Asheim, D., Huybrechs, Wang 2014. Large n asymptotic analysis in Suetin 2009, also D. 2014. Alfredo Dea˜ no Highosc and OPs

Complex Gaussian quadrature Drawbacks: w ( x ) is not a positive weight function, existence of orthogonal polynomials p ω n ( x ) unclear. Location of the quadrature nodes? Expensive calculation (for now...) Advantages: If successful, this approach is optimal and uniform in ω . It is potentially useful for other problems, like integrals with coalescing stationary points or weights of the form w ( x ) = e − V ( z ) , z ∈ C . Related to recent work by Huybrechs, Kuijlaars, Lejon, Mart´ ınez-Finkelshtein, Rakhmanov, Silva... Alfredo Dea˜ no Highosc and OPs

Complex Gaussian quadrature We compute � 1 e i ωx 3 + x d x. − 1 Figure: Errors in log 10 scale with n = 4 (violet), n = 6 (crimson) and n = 8 (black). Alfredo Dea˜ no Highosc and OPs

Where are the zeros? Im � x � 1.0 0.8 0.6 0.4 0.2 1.0 Re � x � � 1.0 � 0.5 0.0 0.5 Figure: Plot of the zeros of p ω n ( x ) with n = 16 as function of ω . Alfredo Dea˜ no Highosc and OPs

Other questions Logarithmic potential interpretation of complex Gaussian quadrature. Computational aspects of numerical steepest descent. Computational aspects of complex Gaussian quadrature. Alfredo Dea˜ no Highosc and OPs

One announcement 2-year Postdoctoral Research Associate Position SMSAS, University of Kent, UK EPSRC First Grant project “Painlev´ e equations: analytical properties and numerical computation” Start date: September 2017 More information: https://www.kent.ac.uk/smsas/ A.Deano-Cabrera@kent.ac.uk Alfredo Dea˜ no Highosc and OPs

The end See you in Canterbury for OPSFA14, July 3-7, 2017 !! https://blogs.kent.ac.uk/opsfa/ Alfredo Dea˜ no Highosc and OPs