EFFICIENT REPRESENTATION AND ACCURATE EVALUATION OF OSCILLATORY INTEGRALS AND FUNCTIONS

. We introduce a new method for functional representation of oscillatory integrals within any user-supplied accuracy. Our approach is based on robust methods for nonlinear approximation of functions via exponentials. The complexity of evaluation of the resulting representations of the oscillatory integrals does not depend or depends only mildly on the size of the parameter responsible for the oscillatory behavior.

1. Introduction.Methods for asymptotic evaluation of oscillatory integrals have a long history (see e.g.[34,10] and references therein).These methods have been widely used to construct asymptotic solutions of partial differential equations (PDEs).Two early examples are [21] as well as [32] where Peter Lax applied an asymptotic approach for solving oscillatory initial value problems.Further development of asymptotic methods led to the theory of pseudo-differential and Fourier integral operators; overall, this work occupied attention of mathematicians for many decades.
While ubiquitous in a variety of applications, computing oscillatory integrals via standard quadrature rules is highly inefficient since the cost of evaluation grows (at best) proportionally to the number of oscillations of the integrand.For example, consider the Fourier-type integral where we assume that the real-valued functions f and g, usually referred to as the amplitude and the phase, are smooth and only mildly oscillatory.The difficulty arises for ω 1 since the integrand becomes highly oscillatory.In order to avoid quadratures, the classical approach to approximate values of I (ω) is to construct its asymptotic expansion with respect to inverse powers of ω.Such asymptotics identifies two main contributions, one from the end points and one from the stationary points (i.e,where g (x) = 0).If g (x) = 0 in [−1, 1], then integration by parts in (1) yields the desired asymptotics.If at an isolated point x * ∈ [−1, 1] one or more derivatives of g vanish, then the asymptotics is obtained using the Taylor expansion of g at x * .The specific powers in the asymptotic expansion depend e.g., on the type of stationary points of g.Asymptotics expansions of this type are also available in higher dimensions, see e.g., [34,10].
More recently, Iserles and Norsett [26,27] developed a Filon-type method for (1) by assuming that the amplitude function f is well approximated by polynomials.The need to solve ordinary differential equations with highly oscillatory forcing terms motivated a combination of asymptotic and numerical approaches taken in [29,15,28,16,17,25].A numerical approach to the evaluation of oscillatory integrals impacts many problems in other applications as may be seen in [19].
In this paper we introduce a new semi-analytic method for the numerical evaluation of oscillatory integrals.Our approach is based on methods for nonlinear approximation of functions via exponentials that yield, for any user-defined accuracy, functional representations of oscillatory integrals.One of the tools we use is the approximation of functions via bandlimited (purely oscillatory) exponentials, an alternative to the traditional approximation by polynomials.Since the integrand of an oscillatory integral has two components, an amplitude and an oscillatory exponential with a large parameter, it is natural and, as we demonstrate, advantageous to approximate the amplitude via exponentials.Indeed, the resulting integrals can be evaluated explicitly yielding a functional representation within any user-selected accuracy.Our construction relies on Gaussian-type quadratures for exponentials described in [4,38].Another tool is the approximation of functions via decaying oscillatory exponentials [5,6,7].We note that the latter methodology also yields near optimal approximations via rational functions.
A combination of these tools allows us to accurately evaluate oscillatory integrals at a cost that does not depend (or depends very mildly) on the size of the parameter ω.Our approach is semi-analytic since it yields a functional approximation, i.e., the result is a (parametrized) function that can be used in further computations.Previously we applied these nonlinear techniques to approximating kernels of operators (see e.g.[7] and references therein) and solving partial differential and integral equations [24,3].We also developed an approach based on nonlinear approximations for applying the oscillatory Rayleigh-Sommerfeld kernel in optics [33].
To demonstrate our approach, we illustrate it on three representative problems.First, in Section 3, we address the accurate computation of integrals of type (1) with complexity either O (log ω) or O (1).Then, in Section 4, we turn to the numerical evaluation of non-traditional oscillatory integrals introduced in [15] (ExpSin integrals), where α > 0, β ∈ R, τ ∈ C, τ = 0.For any user-supplied accuracy, our approach yields a (semi-analytic) representation of this integral such that the number of operations for its evaluation does not depend on ω.It turns out that even the more general integrals in [25] submit to our approach.In Section 5, we present examples of accurate and efficient representation of several (familiar) functions defined in terms of oscillatory integrals and conclude the paper in Section 6.

2.
Preliminaries.Our approach relies on algorithms for representing functions as linear combinations of exponentials with purely imaginary exponents, or decaying oscillatory exponentials, rational functions or Gaussians with complex-valued exponents [4,5,6,23,38,24,3].One of the critical aspects of these algorithms is that they yield representations with near optimal (smallest) number of terms needed to represent a function within a prescribed accuracy.As it becomes clear later in the paper, the very fact that we use exponentials rather than polynomials to approximate amplitudes in oscillatory integrals is an essential advantage of our approach.
2.1.Approximation of functions via exponentials with purely imaginary exponents.While periodic band-limited functions may be expanded into Fourier series, neither the Fourier series nor the Fourier integral may be used efficiently for non-periodic functions on intervals.This motivates considering the linear space of functions with a fixed parameter c, the so-called bandlimit.It is shown in [8] that E c is dense in the space of bandlimited functions, In our approach, given a finite accuracy , we represent functions in E c using a fixed set of exponentials e icθ k x M k=1 , where M = M (c, ).It turns out that by finding quadrature nodes {θ k } M k=1 and weights {w k } M k=1 for exponentials with bandlimit 2c and accuracy 2 , we in fact obtain an approximate basis for E c with accuracy (see [4]).We note that any function on [−1, 1] that is not bandlimited but can be approximated by a bandlimited function with accuracy , can be also represented via the approximate basis for some bandlimit c.
The construction of required quadrature nodes and weights (the so-called generalized Gaussian quadratures for exponentials) can be found in [4] (see also [44] and [38] for different constructions) and may be summarized as Lemma 2.1.For c > 0 and any > 0, there exist nodes where the number of nodes, The nodes and weights maintain the natural symmetry, t l = −t L−l+1 and w l = w L−l+1 .
Remark 1.The construction in [4,38] yields quadratures for band-limited exponentials integrated with a weight w l e ict l x ≤ , where the weight function w is real-valued (in fact, it does not have to be signdefinite, see [38]).If the weight function w is 1 as in Lemma 2.1, then the approach in [4] identifies the nodes of the generalized Gaussian quadratures in (3) as zeros of the Discrete Prolate Spheroidal Wave Functions (DPSWFs) [40], corresponding to small eigenvalues.The size of the eigenvalue determines the accuracy of the quadrature, .In practice, we have tabulated quadratures described in Lemma 2.1 for fixed accuracy, ≈ 10 −15 as well as for ≈ 10 −7 , and organized them by the number of nodes L rather than via the corresponding bandlimit c.
In order to construct an approximate basis for E c , we use Lemma 2.1 to obtain a quadrature for exponentials of bandlimit 2c > 0 and accuracy 2 > 0, yielding M nodes {θ m } M m=1 and weights {w m } M m=1 , so that In [4] these nodes and weights are computed by solving the approximation problem where R m (θ) = M l=1 r lm e icθ l θ , M m=1 r lm e icθmθ l = δ ll are interpolating functions on the nodes {θ m } M m=1 .Then we have The proof of this Lemma as well as a corrected version of the L ∞ estimate found in [4] are presented in the appendix.However, we note that the direct numerical evaluation of the interpolation error indicates that these estimates are somewhat pessimistic.
Remark 2. It may appear that working within standard double precision arithmetic (≈ 16 digits), the accuracy of approximation in (8) is limited to ≈ 10 −8 due to the accuracy requirement in (6).However, by first computing nodes and weights in extended precision ≈ 10 −32 , these nodes and weights can then be used within double precision arithmetic to achieve ≈ 10 −16 in (8) .

2.2.
Prolate spheroidal wave functions.We briefly summarize some of the relevant results in [42,30,31,39,41] (for details see [44,37]).Let us define operators where c > 0 is the bandlimit.The eigenfunctions ψ 0 , ψ 1 , ψ 2 , • • • of Q c coincide with those of F c , and the eigenvalues µ j of Q c are related to the eigenvalues λ j of F c as While all µ j < 1, j = 0, 1, . . ., for large c, the first approximately 2c/π eigenvalues µ j are close to 1.They are followed by a transition region consisting of O(log c) eigenvalues which decay exponentially fast; the rest of the eigenvalues µ j are very close to zero.The functions ψ j are also eigenfunctions of a differential operator [42].In some respects, PSWFs are strikingly similar to orthogonal polynomials, e.g., they are orthonormal and constitute a Chebychev system.
2.3.Numerical construction of interpolating bases for band-limited functions.We briefly review the construction of interpolating functions R m in [4].
Given nodes and weights in (6), we use ( 6) to replace the kernel in (10) and apply (9) to obtain the algebraic eigenvalue problem, Solving ( 12), the approximate PSWFs on [−1, 1] are then defined consistent with (9) as where η j are the eigenvalues and Ψ j (θ l ) the eigenvectors in (12).We then define the interpolating basis for band-limited functions as where M m=1 r lm e icθmθ l = δ ll or, alternatively, Combining last identity with (13), we can write the interpolating function in terms of the approximate PSWFs as, The interpolating functions R k play the same role as the Lagrange interpolating polynomials defined on the Gauss-Legendre nodes.We note that the interpolating functions are also used in the construction of a symplectic, bandlimited collocation implicit Runge-Kutta (BLC-IRK) method for solving ordinary differential equations [9].
In our approach we want to obtain an approximation on [−1, 1] of a bandlimited function f via a linear combination of exponentials e icθ k x M k=1 , However, the direct computation of the coefficients c m is an ill-conditioned problem and we first approximate f using interpolating functions to obtain Since the coefficients r kl in ( 14) are precomputed, this sequence of steps avoids the numerical difficulties of finding the coefficients in ( 17) directly.In this paper we always assume that ( 17) is available for a target accuracy .
2.4.Approximation of functions via decaying oscillatory exponentials.In [5], for a smooth function f (x) and given accuracy > 0, we solve the approximation problem of finding the minimal number of complex coefficients w m and exponents η m (of positive real part) such that For functions singular at x = 0, we formulate (18) on the interval [δ, a], where δ > 0 is a small parameter and replace absolute error by relative error.The function f may be oscillatory, periodic, or non-periodic and we circumvent the constraints of Fourier analysis by optimizing the value of both the exponents and the coefficients, which are now complex-valued.These approximations have significantly fewer terms than Fourier representations or more general constructions like those of the type (17).Our results on exponential approximations have a dual (Fourier) version as approximations by rational functions.To see why, define for any x ≥ 0 and g(x) = g(−x) for negative values of x.The function g(x) is infinitely differentiable everywhere except at x = 0 and its Fourier transform is a real-valued rational function that can be derived analytically, For more details on these algorithms we refer the reader to [5,6,23,24,3].
2.5.Reduction algorithms.Our algorithms seek to find (nonlinear) approximations to a function on a given interval using an optimal (minimal) number of terms for a target accuracy.We found that it is often advantageous to first construct a suboptimal representation of the desired form and then obtain the optimal one via an alternative (reduction) algorithm, which minimizes the number of terms.Typically, an accurate but suboptimal approximation may be relatively easy to obtain e.g., by using a quadrature rule in an integral representation of a function.It is also the case in the context of constructing representations of oscillatory integrals, where reduction algorithms (see [5,7,23]) allow us to first use inefficient but accurate quadratures to provide an initial approximation to then be followed by a reduction step.Importantly, the recently developed reduction algorithm [23] is not only fast but also allows to maintain high relative accuracy.We recently demonstrated the efficiency and accuracy of this algorithm in [24] by solving Burgers' equation with a very small viscosity.
3. Representation of Fourier-type integrals.As an example of the straightforward use of our techniques, we consider a linear phase g(x) = x in (1) and compute assuming that f is well approximated by bandlimited exponentials with bandlimit c, where c ω.As in (17), for a target accuracy , we construct the approximation which immediately gives the explicit approximation Here the number of terms, M , is proportional to the bandlimit c and, therefore, the integral in ( 21) can be efficiently evaluated for any parameter ω at a cost independent of its size.

3.1.
A representative example.This example illustrates our approach not only for a linear phase function g but also in the case of a nonlinear phase as discussed below.We select where b = 120.9513171632071.The choice of parameters n and l is explained in the remark below.This function is displayed in Figure 1 and will be reinterpreted via a change of variables in the next section.We represent f via ( 22) with accuracy 1.5 • 10 −7 , which is consistent with the use of quadratures providing integration accuracy ≈ 10 −15 (see Section 2.3).We then represent the integral I (ω) via (23) and display, in Figure 2, the resulting approximation in the intervals [0, 100] and [10000, 10100].
Remark 3. We construct a bandlimited approximation of the function f on [−1, 1] adequate for the choice n = 4 (n later becomes the order of the zero of the phase function).If we were to choose a larger n in (24), then the denominator would change rapidly near x = −1.In this case we would simply split the interval of integration into appropriate subintervals (controlled by the parameter l in ( 24)) to maintain a reasonable bandlimit and, hence, a reasonable number of terms in the approximation.This is a general strategy in situations where the function to be approximated has very different behaviors across the interval of integration.
Next we turn to the case of a nonlinear phase function g.

3.2.
Oscillatory integrals on an interval without stationary points.We consider (1) for a phase that satisfies g (x) > g 0 > 0 on [a, b].Changing variables, we have g(−1) is moderate, the bandlimit of the integrand's amplitude will not increase significantly.Note that we can always subdivide the interval further so that this ratio remains moderate.
Our example above in Section 3.1 (up to a constant and a phase factor) may be interpreted as a change of variables in the integral where ω = (2 n − 1) /2 n(l+1)+1 ω, so that we reduce the computation of this integral to the case of linear phase.Such change of variables is used, on appropriate subintervals, in the next section.

3.3.
Oscillatory integrals on an interval containing stationary points.The principle of stationary phase [34,10,43] implies that, for large ω, the main contributions to the value of I(ω) in (1) come from either the endpoints of the interval of integration or the stationary points.We say that x * is a stationary point of order n − 1 if the first n − 1 derivatives of g vanish at x * , Iserles and Norsett had shown [26,27] how to build quadratures for I(ω) in (1) by approximating the amplitude f via Hermite polynomial interpolation.Their x k e iωg(x) dx, which are assumed to be known.Also, to avoid computation of derivatives of f , a derivative-free variant is presented in [26].
We demonstrate how to construct a functional representation of such integrals on the canonical example where f is only mildly oscillatory and n ≥ 2 is an integer.The only stationary point of this integral is at x * = 0 and we subdivide the original interval to isolate the stationary point within a sufficiently small interval.We subdivide the interval as follows, so that we approach the stationary point in a hierarchical fashion.The parameter L describing the number of levels of subdivision is chosen later.On all subintervals, except the one about zero, we perform a change of variables in order to use (25).We show below that, since the intervals become smaller when approaching the stationary point x * = 0, the bandlimit of the integrand decreases exponentially fast.Once we reach a sufficiently small bandlimit, we evaluate the integral over −2 −L−1 , 2 −L−1 directly.Hence, by first fixing the desired range of values of ω, the cost of evaluation depends only logarithmically on the maximum size of ω, i.e., it is proportional to the number of levels L in (27).Since the intervals in ( 27) are symmetric about zero, we discuss only those where x > 0. Denoting the change of variables y = 2 n(l+1)+1 x n − (2 n + 1) / (2 n − 1) yields Hence, for any target accuracy, we can always find a value of L such that the contribution of I (n) l (ω) for l > L is negligible.We note that the bandlimit of the exponential e i (2 n −1) 2 n(l+1)+1 ωy in (28) decreases exponentially fast as the parameter l increases.
As in Section 3.1, we approximate and obtain Remark 4. If the power n or the bandlimit of are relatively large, then the number of terms in the approximation of each (ω) may be significant.To reduce the number of terms in our representation of we can split the approximation of f (y) from the one for y . We approximate f in the form (17) and use the results in [7] to efficiently approximate where ρ m and η m are positive parameters and M is of moderate size, even if n is very large.The advantage of this approach relies on our ability, via the reduction algorithms described in Section 2.5, to reduce the overall number of exponentials in the approximation of the product of f (y) and y As a result, the value of I l (ω) can be approximated, with an error bounded by ˜ /2 l+1 , by a linear combination of integrals of the form for some values θ ∈ [−1, 1] and η > 0. Therefore, for each n, l, and target accuracy , there is an ˜ such that the error of approximating I l (ω) is within our target accuracy.
Remark 5.An alternative approximation of the integral in ( 28) can be obtained by constructing quadratures for bandlimited exponentials with a weight function.First, we approximate f in the form (17) (without any additional factor, as above), thus reducing the problem to the computation of integrals of the type where a > 1, α ∈ (0, 1) , and p > 0. Second, this last integral is a particular case of band-limited exponentials integrated with a weight function (see ( 4)) for which accurate quadratures can be obtained.In order to construct these quadratures, it is advantageous (see [4,38]) to obtain the function F (α, p, a) explicitly.To this end, we write substitute in (32) to obtain Using [35, 8.6.4],we arrive at where Γ (α, z) is the incomplete Gamma function.
4. Representations of more complicated oscillatory integrals.

ExpSin integrals (integrals with highly oscillatory periodic components).
Next we develop a simple approach to compute where α > 0, β ∈ R, τ ∈ C, τ = 0. Clearly, the apparent difficulties in computation of this integral arise if ω 1.In this case, the integrand is oscillatory due to the fact that, for large ω, ω (αx + β) covers many periods of the sine function.These type of integrals are considered in [15] (dubbed there ExpSin integrals) in order to develop an accurate method for solving ordinary differential equations with highly oscillatory forcing terms (the interest in such integrals stems e.g., from problems in circuit design).The approach in [15] relies on a combination of constructing the asymptotics of this integral and using quadrature formulas.Our approach for evaluating (33) for any ω ≥ 0 appears to be significantly simpler than that in [15].
In order to estimate the bandlimit of the integrand in the representation of the functions u m , it is enough to estimate the bandlimit of the function h (y) = e τ sin(qy) , for q ∈ (0, π], τ ∈ C, τ = 0, and y ∈ [−1, 1].Using the expansion [1, 9.6.33]with z = −τ and t = ie iqy , we obtain where I n is a modified Bessel function of order n, I n (τ ) = i −n J n (iτ ).Therefore, for accuracy , it is enough to find n 0 > 0 such that yielding q • n 0 as the estimate for the bandlimit.From the asymptotic expansion [1, 9.3.1]for large orders n, we obtain Using Stirling's formula, we conclude that which we use to determine the value of n 0 .
4.1.1.Example.We illustrate the computation of (33) for f (x) = sin (ax) / (ax), with a = 20.42035224833366.For this function we constructed the representation (34) with 20 terms and accuracy 0.2 • 10 −15 .In this example, we set α = 1, β = 7.7 and τ = 2 and check the accuracy of our algorithm using Mathematica T M (with high working precision) to numerically evaluate the integral (33) for several selected values of ω in order to verify the estimate (35).We illustrate the result in Figure 3.
In this example, we evaluated the integral at 10, 000 points in each interval taking about 3.5 seconds (in all intervals) on a laptop (with no attempt to optimize the Fortran 90 code).
Our approach for computing (33) is also applicable to the evaluation of the more general integral
We have with p and q as in (38).The integrals v m and v 0 are only mildly oscillatory and are evaluated using quadratures in [4].We observe that using band-limited exponentials allows us to separate the Dirichlet factor in (39) to combine the contribution of all subintervals except the two small (leftover) subintervals contributions which are evaluated separately.This strategy reduces the problem to the evaluation of three mildly oscillatory integrals.

5.
Representations of oscillatory functions.In this section we discuss additional examples illustrating the efficient representation of highly oscillatory functions.As a first example, consider the Fourier transform of the characteristic function of an ellipse where a and b are positive.In polar coordinates we have If a = b, the function f (ρ, θ) is highly oscillatory not only as a function of ρ 1 (with its bandlimit growing linearly with ρ), but also as a function of θ.Due to the slow decay of f and its rapid oscillations in θ, representations of this function via bases (for example, using curvelets [12,14] or as in [13,18]) are not efficient for large ρ since the number of terms grows quadratically with the bandlimit.On the other hand, using the results in [5,7], we have where M = O(log c) and M = O(log −1 ).In fact, (42) is valid in [0, ∞) once the bandlimit c is chosen to be sufficiently large.We illustrate the error of this approximation in Figure 4.
We note that a representation of f in a basis would require a large number of terms essentially dictated by Nyquist sampling criterion.In contrast our nonlinear approximation (43) circumvents Nyquist constraint and only requires a small number of terms, M .
5.1.One-dimensional oscillatory integral transforms.We now consider the Fourier transform of a radial function f it is easy to see (e.g., by Bochner's theorem [20, pp. 247]), that the univariate function u (ρ) is obtained via the Hankel transform, where J d 2 −1 is the Bessel function of order d 2 − 1 and ρ ≥ 0. We note that if f has singularities, then the decay of u is slow.Writing u as where α = d/2 − 1 and f (t) = (2π) d 2 f (t)t d−1 , we observe that the kernel (ρt) −α J α (ρt) is an oscillatory function.Instead of discretizing (45), we will approximate both, the function f and the kernel by short sums of exponentials.As a consequence, we will obtain a rational representation for the function u(ρ).
First, by an analysis similar to the one in [2, p. 203], we express the kernel function x −α J α (x) as a Laplace type integral, Using the nonlinear algorithms described in [5,7], we can find a contour Γ and the corresponding parametrization γ where the integrand is only mildly oscillatory.Therefore, once the path is selected, it is sufficient to use the trapezoidal rule to discretize the integral to any desired accuracy.Note that the oscillatory behavior of the Bessel function is now encoded in the imaginary part of the exponents τ m in (46).The fact that the trapezoidal rule over the whole real line could be very accurate and efficient follows from the following result.
Let R(h) be the error of using the trapezoidal rule over R with step size h > 0 and shift d ∈ (0, h), A fast decay of the function g implies that, for a given accuracy, only a small number of terms is used in (47).
From [22,Thm,5.13],we have Therefore, when g is analytic and bounded in a strip about the real axis, the trapezoidal rule error decays exponentially fast with the step size h, which is the behavior we have observed for many integrals related to special functions.
Using (46), we rewrite (45) as The number of terms in this expression can be minimized using the results in Section 2.5.Alternatively, using the integral representation in (46), we obtain which could be discretized and optimized via the reduction algorithm of Section 2.5 to yield the desired rational approximation of u.We note that the results for radial functions can be incorporated into a more general construction using the Funk-Hecke formula for the Fourier transform of the

Figure 2 .
Figure 2. Approximation of the integral (21) for the amplitude f in (24) in the intervals [0, 100] (top) and [10000, 10100] (bottom).Real part of f is displayed with dashes, imaginary part with dots and absolute value with a solid line.