# An adaptive Levin method for complicated domains

Shukui Chen$\mbox{}^{\dagger}$, Kirill Serkh$\mbox{}^{\dagger\,\diamond}$,
James Bremer$\mbox{}^{\diamond}$

###### Abstract

In this paper we describe an adaptive Levin method for numerically evaluating
integrals of the form
$\int_{\Omega}f(\mathbf{x})\exp(ig(\mathbf{x}))\,d\Omega$ over general domains
that have been meshed by transfinite elements. On each element, we apply the
multivariate Levin method over adaptively refined sub-elements, until the
integral has been computed to the desired accuracy. Resonance points on the
boundaries of the elements are handled by the application of the univariate
adaptive Levin method. When the domain does not contain stationary points, the
cost of the resulting method is essentially independent of the frequency, even
in the presence of resonance points.

$\mbox{}^{\dagger}$Department of Computer Science and
$\mbox{}^{\diamond}$Department of Mathematics, University of Toronto

##  1 Introduction

We consider multivariate oscillatory integrals of the form

$\int_{\Omega}f(\mathbf{x})\exp(ig(\mathbf{x}))\,d\Omega$ (1)

where $f$ and $g$ are smooth, slowly varying functions, and $g$ is real-
valued. The integration domain $\Omega$ is a bounded region in
$\mathbbm{R}^{2}$ with a piecewise smooth boundary. We do not impose any
additional conditions on $\nabla g$, and allow for stationary points and
resonance points in $\overline{\Omega}$. There exist several approaches for
the numerical evaluation of integrals of the form (1), including Filon-type
methods [1], Levin-type methods [2] and numerical steepest descent [3].

We demonstrate an adaptive Levin method for numerically evaluating the
integral (1) over general domains that have been meshed by transfinite
triangular elements. By computing the integral as the sum of contributions
from the triangular elements, our method can handle arbitrarily complicated
domains, and can be run in parallel over the elements. Over each triangular
element, we apply the multivariate Levin method on adaptively refined sub-
elements of the triangular element, until the integral has been computed to a
desired accuracy. Resonance points on the boundaries of the elements are
handled by the application of the univariate adaptive Levin method.

##  2 The Levin PDE

The multivariate Levin method operates by solving the partial differential
equation

$\mathcal{L}[\mathbf{p}]:=\nabla\cdot\mathbf{p}+i\nabla g\cdot\mathbf{p}=f$
(2)

in order to find a vector field $\mathbf{p}$ such that

$\nabla\cdot(\mathbf{p}(\mathbf{x})\exp(ig(\mathbf{x})))=f(\mathbf{x})\exp(ig(\mathbf{x})).$
(3)

By applying the divergence theorem, the integral (1) is transformed into

$\int_{\Omega}f(\mathbf{x})\exp(ig(\mathbf{x}))\,d\Omega=\int_{\Omega}\mathcal{L}[\mathbf{p}](\mathbf{x})\exp(ig(\mathbf{x}))\,d\Omega=\int_{\Gamma}\mathbf{n}(\mathbf{x})\cdot\mathbf{p}(\mathbf{x})\exp(ig(\mathbf{x}))\,d\Gamma,$
(4)

where $\Gamma$ is the boundary of the integration domain. The univariate
integral over the boundary $\Gamma$ can then be evaluated using the univariate
adaptive Levin method.

It was long believed that this approach suffers from “low-frequency
breakdown,” meaning that the accuracy of the calculated value of the integral
deteriorates when the integrand is slowly varying, as in the case of a
stationary point of $g$. However, it was proven in the univariate case that if
a Chebyshev spectral method is used to discretize the differential equation
and the resulting linear system is solved via a truncated singular value
decomposition, then no low-frequency breakdown occurs [4]. The proof
guarantees that in the univariate case, high accuracy can be obtained
regardless of the magnitude of $g^{\prime}$ and whether or not it has zeros.
The same argument can be readily generalized to the multivariate case to show
that stationary points of $g$ need not be a restriction when applying the
multivariate Levin method adaptively.

##  3 Numerical Discretization

To solve the Levin PDE numerically, we seek a vector field $\mathbf{p}$ of the
form $\mathbf{p}(\mathbf{x})=\sum_{j=0}^{n}p_{j}\varphi_{j}(\mathbf{x})$ where
$\\{\varphi_{0},\dots,\varphi_{n}\\}$ is the set of all monomials of degree
less than or equal to $k$, so that $n+1=\text{dim}(\mathcal{P}_{k})$. Using
the collocation method, we solve the system

$\sum_{j=0}^{n}\mathcal{L}[\varphi_{j}](\mathbf{x}_{k})p_{j}=f(\mathbf{x}_{k}),\quad
k=0,\dots,m,$ (5)

for a set of collocation nodes $\\{\mathbf{x}_{0},\dots,\mathbf{x}_{m}\\}$,
where $m+1=\text{dim}(\mathcal{P}_{\ell})$, see [5]. The resulting linear
system is then solved using a truncated singular value decomposition.
Consistent with our analysis in the univariate case, the discretized PDE can
be solved with better accuracy if the collocation nodes have slightly higher
order than the basis, so that $\ell>k$.

By incorporating phase function methods for ordinary differential equations,
the multivariate Levin method can be applied to a large class of oscillatory
integrals involving special functions, including products of such functions
and compositions of such functions with slowly-varying functions, without the
need for any symbolic or arbitrary precision computations [4].

##  4 Resonance Points

A point of resonance is any point $\mathbf{x}\in\Gamma$ such that the vector
$\nabla g(\mathbf{x})$ is orthogonal to the boundary [6]. The resonance points
represent stationary points of the oscillator in the boundary integral (4),
which means that the oscillatory boundary integral cannot be efficiently
evaluated using the non-adaptive univariate Levin method of fixed order
without sampling at the stationary points. A nonresonance condition is usually
imposed to require the vector $\nabla g(\mathbf{x})$, $\mathbf{x}\in\Gamma$,
to be nowhere orthogonal to the boundary.

We note that the vector field $\mathbf{p}$ can often be found to high
accuracy, in the sense that the error
$\|\mathcal{L}[\mathbf{p}]-f\|_{L^{\infty}(\Omega)}$ is small regardless of
whether or not resonance points are present. This suggests that the integral
(1) can be accurately approximated by the boundary integrals (4). Therefore,
we apply the univariate adaptive Levin method to efficiently evaluate the
boundary integrals, since it can handle the case of stationary points
automatically. As a consequence, we can evaluate of the integral (1)
efficiently by applying the Levin method adaptively and recursively,
regardless of the magnitude of $\nabla g$ and whether stationary points or
resonance points are found in $\overline{\Omega}$.

##  5 Numerical experiment

We demonstrate the efficiency and accuracy of our method by computing the
integral

$\int_{\Omega}\nabla\cdot(G(\mathbf{x},\mathbf{x}^{\prime})\nabla
u(\mathbf{x}^{\prime}))\,d\mathbf{x}^{\prime},\quad\mathbf{x}=\mathbf{0},$ (6)

where
$G(\mathbf{x},\mathbf{x}^{\prime})=\frac{i}{4}H_{0}^{(1)}(\omega\|\mathbf{x}-\mathbf{x}^{\prime}\|)$
is the Green’s function of the 2D Helmholtz equation and
$u(\mathbf{x}^{\prime})=\|\mathbf{x}^{\prime}\|_{2}^{2}$ is chosen to be a
smooth function. We use the non-oscillatory phase function for Bessel’s
equation to transform the integral (6) into the form (1).

The integration domain $\Omega$ (Figure 1) was constructed and tessellated
using Gmsh and its built-in CAD kernel [7]. After tessellation, we constructed
piecewise Chebyshev expansions to represent the boundary as curved segments,
and built the geometric mapping using transfinite interpolation [8].

Figure 1: Real part of the integrand of (6) with $\omega=100$.

We point out that, due to the geometry, every point on the curved boundary in
the bottom left is a resonance point. Moreover, due to the tessellation,
resonance points are also found on some edges of the interior triangles.
However, our method was nonetheless able to evaluate the integral (6) in a
nearly frequency-independent manner.

For estimating the absolute error obtained by our method, we transformed the
integral into a boundary integral by noting that the integrand is the
divergence of a vector field. We then evaluated the boundary integral
adaptively using a 30-point Gauss-Legendre quadrature rule.

##  6 Numerical experiment

We present in Figure 2 the results of the numerical experiments computing the
integral (6), in order to illustrate the properties of the multivariate
adaptive Levin method. The code for the experiment was written in Fortran and
compiled with version 13.2.0 of the GNU Fortran compiler. The experiments were
performed on a laptop computer equipped with an Intel i9-12900HK processor and
32GB of memory. We used OpenMP to parallize over the triangular elements.

Figure 2: The results of the numerical experiments computing the integral (6).
The top plot gives the error in the value of the integral computed via the
multivariate adaptive Levin method as a function of $\omega$. The bottom plot
gives the running time in milliseconds as a function of $\omega$.

##  7 References

## References

  * [1] J. Gao and A. Iserles, “On an extended Filon method for highly oscillatory integrals over a simplex,” Mathematics of Computation, vol. 92, pp. 867–893, Mar. 2023.
  * [2] S. Olver, “On the Quadrature of Multivariate Highly Oscillatory Integrals Over Non-polytope Domains,” Numerische Mathematik, vol. 103, pp. 643–665, June 2006.
  * [3] D. Huybrechs and S. Vandewalle, “The Construction of cubature rules for multivariate highly oscillatory integrals,” Mathematics of Computation, vol. 76, pp. 1955–1980, Oct. 2007.
  * [4] S. Chen, K. Serkh, and J. Bremer, “On the adaptive Levin method,” Jan. 2024\.  arXiv:2211.13400 [cs, math].
  * [5] B. Vioreanu and V. Rokhlin, “Spectra of Multiplication Operators as a Numerical Tool,” SIAM Journal on Scientific Computing, vol. 36, pp. A267–A288, Jan. 2014.
  * [6] A. Iserles and S. P. Nørsett, “On the computation of highly oscillatory multivariate integrals with stationary points,” BIT Numerical Mathematics, vol. 46, pp. 549–566, Sept. 2006.
  * [7] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009.
  * [8] A. Perronnet, “Interpolation transfinie sur le triangle, le tétraèdre et le pentaèdre. Application à la création de maillages et à la condition de Dirichlet,” Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, vol. 326, pp. 117–122, Jan. 1998.