# Numerical evaluation of highly oscillatory integral

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?

• Is the function even? Jun 4 '19 at 1:40
• Yes, it is even Jun 4 '19 at 5:59
• Have you tried turning your integral into an ODE? Jun 4 '19 at 11:49
• No, differentiation w.r.t. $x$ and then solve the differential equation numerically. Jun 4 '19 at 15:35
• Your first plot seems to show a different function than your integrand. Namely, it seems to have $\lambda$ replaced with $\lambda x$. I.e. the plot is of the function $x\mapsto(\cos(\lambda x\cos x)\operatorname{sinc} x)$ Aug 4 '19 at 9:12

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.

• Thanks, this is a very good idea! Jun 4 '19 at 21:43

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+f)/2+(f+f)/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.

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.

• Nice! Are the instructors of the course your real-life professors? Their course is fantastic, though very tough and fast-paced Jun 5 '19 at 14:10
• @doetoe yes, I am a student of Konstantin. He shared with me a link to your question. Jun 5 '19 at 14:32

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
oct       = 0.943661538066060288748574485680878906503533004997613278231689169604876
asymptotic= 0.944029734



Here's the code:

#include <iostream>
#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;
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$$: • Thanks, this is really nice! I didn't get it to work yet, my boost installation is not compatible, but I'm downloading the latest version now. Jun 4 '19 at 21:37
• Just to be sure: in 23 you have cos(lambda*cos(x))/x, without the sin(x) factor from the integrand. Is it ooura_fourier_sin that assumes this factor sin(x) to multiply the integrand passed to it? Jun 4 '19 at 21:55
• I got it working. It and its dependencies seem to be all header only, so I didn't even have to install or compile (except for the executable). I hope it will get included in boost! Jun 4 '19 at 22:38
• @doetoe: Yes, the $\sin(x)$ factor is implicit and incorporated into the quadrature weights. I'm considering getting it speed up by refactoring it to cache the nodes and weights; we'll see! Jun 4 '19 at 23:06
• @doetoe: It's merged into Boost 1.71. The API is a bit different than this answer gives. Aug 4 '19 at 13:29