Asymptotics for Orthogonal Polynomials and High-frequency Scattering Problems

The goal of this thesis is to exploit asymptotic behaviour in high-degree orthogonal polynomials and high-frequency acoustic scattering problems to obtain %[Dave] problems. This leads to
a lower computational cost. The code is made publicly available and validates this goal as well as the accuracy of the results of this thesis.

We are interested in the higher-order asymptotic behavior of orthogonal polynomials of Jacobi-, Laguerre- and Hermite-type as their degree n goes to ∞. These are orthogonal with respect to the weight functions w(x)  = (1-x)^α (1+x)^β h(x) for x ∈ [-1,1],  w(x) = x^α exp[-Q(x)] for x ∈ [0,\infty), and w(x) = exp[-Q(x^2)] for x ∈ (-\infty, \infty), respectively, with α and β > -1. The functions h(x) and log Q(x) are real analytic and strictly positive on the interval of orthogonality. This information is implicitly available in two papers, where the authors use the Riemann-Hilbert formulation and the Deift-Zhou non-linear steepest descent method. We show that the computation of higher-order terms can be simplified analytically in a jump relation, leading to their efficient, fully explicit and automatic construction. For Laguerre- and Hermite-type polynomials, we need the derivation of an asymptotic expansion of a zero of a general polynomial with respect to its constant coefficient to obtain more explicit expressions when Q(x) is a polynomial. The resulting asymptotic expansions in four different regions of the complex plane are implemented both symbolically and numerically, along with heuristics to determine these regions and the number of terms to use for a certain accuracy.

The main advantage of these expansions is that the accuracy improves as the degree of the polynomials rises, at a computational cost that is practically independent of the degree or even decreasing due to requiring fewer terms. They allow more accurate approximations for modest degrees n than the leading order term (or some higher order terms) that already existed in the literature. In contrast, the typical use of the recurrence relation for orthogonal polynomials in computations leads to a cost that is at least linear in the degree and a decreasing accuracy.

A smooth integral can be approximated with a weighted sum of n integrand evaluations which is stable and optimal for polynomial integrands. The nodes of these Gaussian quadrature rules are zeros of the corresponding orthogonal polynomials corresponding to the weight function in the integral. Those can be computed via a Newton method on the recurrence relation in O(n^2), or on the asymptotic expansions mentioned before. In order to reduce the constant for the latter O(n) algorithm, to reduce code length and to solve possible convergence as well as other numerical issues, we set up explicit asymptotic expansions of the nodes and weights of Gauss-Laguerre and -Hermite quadrature rules with a high order in the relevant regions.

In another part of our research, Boundary Element Methods are suited to solve wave propagation and scattering problems in acoustics, as they only discretise the boundary of the acoustic domain to obtain linearised sound waves in the frequency domain. They lead to a discretisation matrix that is typically dense and large: its size and condition number grow with increasing frequency. Yet, high frequency scattering problems are intrinsically local in nature, which is well represented by highly localized rays bouncing around. Asymptotic methods can be used to reduce the size of the linear system, even making it frequency independent, by explicitly extracting the oscillatory properties from the solution using ray tracing or analogous techniques. However, ray tracing becomes expensive or even intractable in the presence of (multiple) scattering obstacles with complicated geometries. Instead, we start from the same discretisation that constructs the fully resolved, large and dense matrix, and achieve asymptotic compression by explicitly localizing the Green function.

This results in a large but sparse matrix, with a faster associated matrix-vector product and, as numerical experiments indicate, a much improved condition number. Although an appropriate localisation of the Green function also depends on asymptotic information unavailable for general geometries, we can construct it adaptively in a frequency sweep from small to large frequencies in a way that automatically takes into account a general incident wave. We approximately maintain the discretisation error and do not use a priori knowledge on the solution.

We show that the approach is robust with respect to non-convex, multiple and even near-trapping domains, though the compression rate is clearly lower in the latter case. It can be applied to general incident waves, boundary conditions, dimensions and integral formulations. Furthermore, in spite of its asymptotic nature, the method is robust with respect to low-order discretisations such as piecewise constants, linears or cubics, commonly used in applications. On the other hand, we do not decrease the total number of degrees of freedom compared to a conventional, classical discretisation. The combination of the sparsifying modification of the Green function with other accelerating schemes, such as the Fast Multipole Method, appears possible in principle. A visibility approach was shown to be inferior, we have extended the principle of Fermat, and regions that are important for insight into the high-frequency behaviour of the problem are adaptively provided.

Finally, we focus on multiple scattering obstacles, where the wave pattern becomes very complicated, but wavenumber-independent simulation schemes have been proposed in the literature based on ray tracing. In such schemes, one can note that the phases of the corresponding densities on each of the obstacles converge to an equilibrium after a few iterations. The equilibrium is given by periodic orbits between a subset of the obstacles, on which the existing analysis of ray tracing schemes focuses. We approximate the phase of such a limiting mode representing a full cycle of reflections with a Taylor series. We exploit symmetry in the case of two circular scatterers, but also provide an explicit algorithm for an arbitrary number of general 2D obstacles. The coefficients, as well as the time to compute them, are independent of the wavenumber and the incident wave. We sketch the principle to obtain a series expansion of the smoothly varying part of this mode using the method of steepest descent. The results may be used to accelerate ray tracing schemes after a small number of initial iterations.