Numerical evaluation of highly oscillatory integral

Question

In this advanced course on applications of complex function theory at one point in an exercise the highly oscillatory integral

$$I(\lambda)=\int_{-\infty}^{\infty} \cos (\lambda \cos x) \frac{\sin x}{x} d x$$

has to be approximated for large values of $\lambda$ using the saddle point method in the complex plane.

Due to its highly oscillatory nature, this integral is very hard to evaluate using most other methods. These are two fragments of the graph of the integrand for $\lambda = 10$ at different scales:

A leading order asymptotic approximation is

$$I_{1}(\lambda) = \cos \left(\lambda-\frac{1}{4} \pi\right) \sqrt{\frac{2 \pi}{\lambda}}$$

and a further (much smaller) refinement adds the term

$$I_2(\lambda)=\frac{1}{8} \sin \left(\lambda-\frac{1}{4} \pi\right) \sqrt{\frac{2 \pi}{\lambda^{3}}}$$

A graph of the approximated values as a function of $\lambda$ looks as follows:

Now comes my question: to visually see how good the approximation is, I'd like to compare it to the "real value" of the integral, or more precisely to a good approximation to the same integral using an independent algorithm. Due to the smallness of the subleading correction, I would expect this to be real close.

I tried to evaluate the integral for some $\lambda$ using other algorithms, but with very little success: Mathematica and Matlab using the default numerical integrator don't manage to produce a meaningful value (and report this explicitly), mpmath using both the doubly exponential $\tanh(\sinh)$ substitution and the Gauss-Legendre method produces very noisy results, though it does have a slight tendency to oscillate around the values that the saddle point method gives, as this graph may show:

Finally I tried my luck with a Monte-Carlo integrator using importance sample that I implemented, but I didn't manage to get any stable results either.

Does anyone have an idea of how this integral could be independently evaluated for any fixed value of $\lambda > 1$ or so?

Answer 1

Use Plancherel's theorem to evaluate this integral.

The basic idea is that for two functions $f,g$,

$$ I=\int_{-\infty}^{\infty} f(x) g^*(x)dx = \int_{-\infty}^{\infty} F(k) G^*(k) dk $$

where $F,G$ are the Fourier transforms of $f,g$. Your functions both have relatively small support in the spectral domain. Here, $\sin x / x \rightarrow \text{rect}(k)$ and $\cos(\lambda\cos x)$ should have an analytical Fourier transform (or series), like the Jacobi-Anger expansion. You can truncate the infinite series at about $\lambda$ terms due to the super-exponential decay of the Bessel function $|J_n(x)|$ for $n>|x|$. Hope this helps.

Edit: Actually, you should use the Fourier series representations here instead of the transforms. The transform path leads to deriving the asymptotic representation you already have (turns out this is just $\pi J_0(\lambda)$). Plancherel's theorem above also works for Fourier series with an integration domain of $[0,2\pi]$ on the last integral.

Answer 2

The key to the evaluation of oscillatory integrals is to truncate integral at the right point. For this example you need to choose upper limit of the form $$ \pi\mathbb{N}+\frac{\pi}{2} $$ Before explaining why it should work, let me at first show that it actually produces good results.

Asymptotics

It is easy to guess that asymptotic series has the form $$ I(\lambda) \sim \sqrt{\frac{2\pi}{\lambda}}\left[\cos\left(\lambda-\frac{\pi}{4}\right) + c_1 \frac{\sin\left(\lambda-\frac{\pi}{4}\right)}{\lambda} + c_2 \frac{\cos\left(\lambda-\frac{\pi}{4}\right)}{\lambda^2} + c_3 \frac{\sin\left(\lambda-\frac{\pi}{4}\right)}{\lambda^3} + \dots \right] $$ In order to check numerically that $c_1 = \frac{1}{8}$ it is enough to plot difference between an integral and leading asymptotic expression.

int := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x, 0, 20.5*Pi}]; 
Plot[{l*(Sqrt[2*l/Pi]*int - Cos[l-Pi/4]), Sin[l-Pi/4]/8}, {l, Pi/4, 20}]

As an output you get quite a nice sine which coincides with the one you derived above.

If you want to find following coefficients, a little bit more sophisticated piece of code if needed. The idea of the code below is to take several high-lying upper limit values and to "average" their results.

J[l_?NumericQ] := Block[{n=500},
  f[k_] := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x,0,(n+k)*Pi+Pi/2},
  Method->{"DoubleExponential"}, AccuracyGoal->14, MaxRecursion->100];
  1/2*((f[0]+f[1])/2+(f[1]+f[2])/2)
]
t = Table[{l, l^2*(Sqrt[2*l/Pi]*J[l] - Cos[l-Pi/4] - 1/8*Sin[l-Pi/4]/l)}, 
    {l, 4*Pi+Pi/4, 12*Pi+Pi/4, Pi/36}];
Fit[t, Table[Cos[l-Pi/4+Pi/2*n]/l^n, {n, 0, 10}], l]

That produces the following answer. $$ \boxed{ \boxed{ c_2 = -\frac{9}{128}, \qquad c_3 = -\frac{75}{1024}, \qquad c_4 = \frac{3675}{32768}, \quad\dots }} $$

Explanation

Simple example

For the illustration I am going to use more simple example of the sine--integral $$ S(x) = \int_{0}^x \frac{\sin(y)}{y} dy. $$ Let me imagine that I am interested in value $S(\infty)=\frac{\pi}{2}$, but I do not know it.

You see that $S(x)$ oscillates around its limiting value similarly to how partial sums of alternating in sign series oscillate with upper cut-off. $$ S_N = \sum_{n=1}^N \frac{(-1)^n}{n}. $$ When you want to estimate such sum, according to Euler series acceleration method you should take $$ S \approx S_N + \frac{1}{2}\frac{(-1)^{N+1}}{N+1}. $$ Or in terms of sine--integral function, you should integrate it up to the point between maximum and minimum of oscillations. As it is clearly seen from the plot such point are approximately given by $$ S(x) \approx \int_0^{\pi\mathbb{N}+\frac{\pi}{2}} \frac{\sin x}{x}dx $$ for large values of argument. More generally, such point is the one where $\max |S'(x)|$ occurs.

Your problem

Coming back to integral from Konstantin and Yaroslav's course you can see that it behaves exactly the same way as sine--integral as a function of upper limit. That means you only need to calculate values $$ I_{x_0}(\lambda) = 2\int_{0}^{x_0}\cos\left(\lambda\cos(x)\right)\operatorname{sinc}(x)dx $$ with $x_0 = \pi\mathbb{N}+\frac{\pi}{2}$. Below is the plot of several such values with $\lambda = 12\pi$.

tab = Table[{x0, 2*NIntegrate[Cos[12*Pi*Cos[x]]*Sinc[x], {x, 0, x0}, 
     Method->{"DoubleExponential"}, AccuracyGoal->12, MaxRecursion->100]},
    {x0, 10*Pi+Pi/2, 30*Pi+Pi/2, Pi}];
tab1 = Table[(tab[[i]] + tab[[i+1]])/2, {i,1,Length[tab]-1}];
ListPlot[{tab, tab1}]

Here you can see the result of another acceleration method. I rearrange partial sums in the following manner $$ S_N' = \frac{1}{2}(S_N+S_{N+1}) $$ and obtain new sequence $S_N'$ which converges much faster. That trick also happens to be useful if you want to evaluate integral with high precision.

Answer 3

Ooura's method for Fourier sine integrals works here, see:

Ooura, Takuya, and Masatake Mori, A robust double exponential formula for Fourier-type integrals. Journal of computational and applied mathematics 112.1-2 (1999): 229-241.

I wrote an implementation of this algorithm but never put in the work to get it fast (by, say caching nodes/weights), but nonetheless, I'm getting consistent results at everything beyond float precision:

float     = 0.0154244
double    = 0.943661538060268
long dbl  = 0.943661538066058702
quad      = 0.943661538066060288748574485677942
oct       = 0.943661538066060288748574485680878906503533004997613278231689169604876
asymptotic= 0.944029734

Here's the code:

#include <iostream>
#include <boost/math/quadrature/ooura_fourier_integrals.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

template<class Real>
Real asymptotic(Real lambda) {
    using std::sin;
    using std::cos;
    using boost::math::constants::pi;
    Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
    Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
    return I1 + I2;
}

template<class Real>
Real osc(Real lambda) {
    using std::cos;
    using boost::math::quadrature::ooura_fourier_sin;
    auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
    Real omega = 1;
    Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega);
    return 2*Is;
}

template<class Real>
void print(Real val) {
   std::cout << std::defaultfloat;
   std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
   std::cout <<  val <<  " = " << std::hexfloat << val;
}

int main() {
    using boost::multiprecision::float128;
    float128  s = 7;
    std::cout << "Asymptotic = " << std::setprecision(std::numeric_limits<float128>::digits10) << asymptotic(s) << "\n";
    std::cout << "float precision = ";
    print(osc(7.0f));
    std::cout << "\n";
    std::cout << "double precision= ";
    print(osc(7.0));
    std::cout << "\n";
    std::cout << "long double     = ";
    print(osc(7.0L));
    std::cout << "\n";
    print(osc(s));

    print(osc(boost::multiprecision::cpp_bin_float_oct(7)));
    std::cout << "\n";
}

You can't really see the difference between the quadrature and the asymptotic because they lie right on top of each other except as $\lambda\to 0$: