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

Quadrature of highly oscillatory integrals: the role

• f (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 1D integral has the form I[f] = b

a

f(x)w(x)dx ≈

n

• k=1

bkf(xk) = Q[f], where w(x) is a positive weight function, with finite moments xk are the quadrature nodes bk 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 xk = a + kh, k = 0, 1, . . . , n, h = b − a 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 ∈ Pm, for m as large as possible. Here Pm = {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

• f exactness equal to d = n − 1 + k if and only if

1

The quadrature is interpolatory

2

The nodal polynomial pn(x) =

n

• k=1

(x − xk) satisfies b

a

pn(x)q(x)w(x)dx = 0, q ∈ Pk−1

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 pn(x): b

a

pn(x)xkw(x)dx = 0, k = 0, 1, 2, . . . , n − 1, 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 I[f] = b

a

f(x)eiωxdx, where ω ≫ 1. Other examples are I[f] = b

a

f(x) sin(ωx) cos(ωx)

• dx,

I[f] = b

a

f(x)H(1,2)

ν

(ωx)dx,

• r more in general

I[f] = b

a

f(x)eiωg(x)dx, 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 I1[f] = 1

−1

1 2 + xeiωxdx, I2[f] = 1

−1

1 2 + xeiωx2dx.

Figure: Errors in log10 scale in approximating I1[f] (left) and I2[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 I[f] = 1

−1

f(x)eiωxdx = −

s

• k=0

1 (−iω)k+1

• f(k)(1)eiω − f(k)(−1)e−iω

+ O(ω−s−1) 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 2 3 4 5 6 7 8 9 10 −1 −0.5 0.5 1 1 2 3 4 5 6 7 8 9 10 −1 −0.5 0.5 1

Figure: Plot of cos(10 − x) sin ωx, ω = 10 (top) and ω = 1000 (bottom).

Alfredo Dea˜ no Highosc and OPs

Stationary points

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.5 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.5 0.5 1

Figure: Plot of cos(100 x2) (top) and cos(100 x3) (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 2s + 1 that satisfies p(i)(−1) = f(i)(−1), p(i)(1) = f(i)(1), i = 0, 1, . . . , s, and then Q[f] = 1

−1

p(x)eiωxdx. We try the same example as before: I1[f] = 1

−1

1 2 + xeiωxdx, I2[f] = 1

−1

1 2 + xeiωx2dx.

Alfredo Dea˜ no Highosc and OPs

Other methods II

Figure: Errors in log10 scale in approximating I1[f] (left) and I2[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)

n

(x). Filon–Clenshaw–Curtis, with zeros of Chebyshev polynomials P

( 1 2 , 1 2 ) n

(x) (Dom´ ınguez, Graham, Smyshlyaev 2011). Adaptive Filon.

Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent

A different alternative is the following: take I[f] = 1

−1

f(x)eiωg(x)dx and deform [−1, 1] following paths of steepest descent of g(x) in C. For example, for g(x) = x:

Re z Im z −1 + iP 1 + iP 1 −1 iP

C

Figure: Steepest descent paths for the oscillatory integral 1

−1 f(x)eiωxdx.

Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent

We choose P > 0 and we write I[f] = 1

−1

f(x)eiωxdx = i P f(−1 + ip)eiω(−1+ip)dp + 1

−1

f(x + iP)eiω(x+iP)dx − i P f(1 + ip)eiω(1+ip)dp

Alfredo Dea˜ no Highosc and OPs

Numerical steepest descent

We choose P > 0 and we write I[f] = 1

−1

f(x)eiωxdx = ie−iω P f(−1 + ip)e−ωpdp + e−ωP 1

−1

f(x + iP)eiωxdx − ieiω P f(1 + ip)e−ωpdp. 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 I[f] = 1

−1

f(x)eiωxdx, and apply Gaussian quadrature directly, with weight w(x) = eiω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

−1

eiωx 3 + xdx.

Figure: Errors in log10 scale with n = 4 (violet), n = 6 (crimson) and n = 8 (black).

Alfredo Dea˜ no Highosc and OPs

Where are the zeros?

1.0 0.5 0.0 0.5 1.0 Rex 0.2 0.4 0.6 0.8 1.0 Imx

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