Another interesting integral related to the Omega constant is the following $$\int^\infty_0 \frac{1 + 2\cos x + x \sin x}{1 + 2x \sin x + x^2} dx = \frac{\pi}{1 + \Omega}.$$ Here $\Omega = {\rm W}_0(1) = 0.56714329\ldots$ is the Omega constant while ${\rm W}_0(x)$ is the principal branch of the Lambert W function. A similar interesting integral can be found here.

When attempting to check the above answer numerically, Mathematica 9.0 gives $2.0101 \pm 0.0005$ (with a complaint about failure to converge to prescribed accuracy after 9 recursive bisections...) compared to the exact answer of $2.004662032\ldots$

Plotting the integrand two observations can be made: (i) it tends to $\infty$ as $x \rightarrow 0^+$, and (ii) for $x$ greater than about 3 it is highly oscillatory with the amplitude of the oscillations tending to zero as $x \rightarrow \infty$.

Given this, what would be the best way to evaluate the integral numerically if an accuracy of say $10^{-9}$ is needed? This accuracy should be achieved in a "reasonable" amount of time.

  • 1
    $\begingroup$ Mathematica can be coaxed into producing a numerical result with an error of $\approx 2.5 \times 10^{-8}$ with the command NIntegrate[(1+2 Cos[x]+x Sin[x])/(1+2 x Sin[x]+x^2),{x,0,3000},WorkingPrecision->15,MaxRecursion->1000,MinRecursion->500,MaxPoints->40000]+NIntegrate[Sin[x]/x,{x,3000,\[Infinity]},WorkingPrecision->20,MaxRecursion->300,MinRecursion->100]. Some of these parameters are probably overkill. $\endgroup$ Mar 5, 2014 at 21:54

1 Answer 1


Denote $$F(x) = \frac{1 + 2\cos x + x \sin x}{1 + 2x \sin x + x^2}.$$ Our aim is to compute $$I = \int_0^\infty F(x) ~ dx.$$ Using numerical quadrature specifically designed for oscillatory functions – quadosc, part of mpmath, in Python 2.7 – I get 50 decimal digits of accuracy in ~34 seconds:

>>> from mpmath import *
>>> import time
>>> mp.dps = 50
>>> mp.pretty = True
>>> F = lambda x: (1+2*cos(x)+x*sin(x))/(1+2*x*sin(x)+x*x)
>>> t=time.time(); I=quadosc(F, [0,inf], period=2*pi); elapsed=time.time()-t;
>>> I
>>> elapsed

Compare with WolframAlpha, or with mpmath:

>>> +pi/(1 + lambertw(1, 0))

What's the underlying math here? In Section 1, I'll write my understanding of the quadosc function, as it applies to a function $$f(x) = g(x) \cos(\omega x + \phi),$$ where $g(x)$ is “slowly decaying.” In Section 2, I'll outline some of the important ideas / give a non-rigorous explanation of how quadosc works for $F(x)$—the “tl;dr” version is that $F(x)$ decays asymptotically and its zeros “asymptotically agree” with $\sin x$.

1. quadosc for product of slowly-decaying function and sine wave

Suppose $f(x) = g(x) \cos(\omega x + \phi)$ for some slowly-decaying function $g(x)$. Without loss of generality, suppose $\phi=0$. In this case, quadosc will break the integral $$I = \int_0^\infty f(x) ~ dx$$ into half-periods: $$I = \sum_{k = 0}^\infty \int_{k(2\pi/\omega)/2}^{(k+1)(2\pi/\omega)/2} f(x) ~ dx. \tag{*}$$ In other words, it will compute the integral on $[0,\pi/\omega]$, $[\pi/\omega,2\pi/\omega]$, $[2\pi/\omega,3\pi/\omega]$, $[3\pi/\omega,4\pi/\omega]$, etc. (From the source code, it looks like quadosc uses Gauss–Legendre quadrature on each subinterval.)

The resulting series (*) is an alternating series! Therefore, nsum – the mpmath function for numerically approximating infinite series – is called. According to its documentation, by default, nsum attemps to compute the sum two different ways – Richardson extrapolation, or Shanks transfomation, which are both sequence/series acceleration methods.

Presumably both of these components – numerical integration on each subinterval, and numerical series summation – are done until the desired digits of precision (mp.dps = 50, where dps = decimal digits of precision) are satisfied.

2. quadosc for $F(x)$

Why does quadosc work for $\int_0^\infty F(x) ~ dx$? Because $F(x)$ asymptotically decays, and because its zeros “asymptotically agree” with $\sin(x)$, in the sense that “the positive roots $(x_n)_{n=1}^\infty$ of $F(x)$ (where $x_n < x_{n+1}$) satisfy $\lim_{n\to\infty} (n\pi - x_n) = 0$.

How? I'll make two claims, and outline the important ideas that (I predict) could be used to derive a rigorous proof. (1) The roots of $1+2\cos x + x\sin x = 0$ “asymptotically agree” with the roots of $1 + x\sin x = 0$. (2) The roots of $1+ x \sin x = 0$ “asymptotically agree” with the roots of $\sin x = 0$

Important ideas for (1): Invoking the asymptotic series of $\sqrt{4+x^2}$ and $\arctan(x/2)$ – computed with Wolfram Alpha [1], [2] – we have \begin{align*} \sqrt{4+x^2} &\sim x + \frac2x - \frac{2}{x^3} + \frac{4}{x^5} + O(x^{-7}) \\ \arctan\frac x2 &\sim \frac\pi2 - \frac2x + \frac{8}{3x^3} - \frac{32}{5x^5} + O(x^{-7}) \end{align*} Therefore, the numerator of $F(x)$ – which determines the roots of $F(x)$ – satisfies \begin{align*} 1+2\cos(x) + x\sin(x) &= 1 + \sqrt{4+x^2}\cos\left(x - \arctan\frac{x}{2}\right) \\ &\sim 1 + \left(x + \frac2x + O(x^{-3})\right) \cos\left( x - \frac\pi2 + \frac2x + O(x^{-3}) \right). \end{align*} Dropping the $x^{-1}$ and higher-order terms, and recalling $\cos(x - \frac\pi2) = \sin x$, $$ 1+2\cos(x) + x\sin(x) \sim 1 + x\sin(x).$$ As $x\to\infty$, the difference between $1+2\cos(x) + x\sin(x)$ and $1 + x\sin(x)$ will vanish, which implies that their roots will asymptotically agree, as claimed.

Important ideas for (2): Note that $$1 + x\sin x = 0 \iff \frac1x + \sin x = 0 \quad (\text{assuming } x\not=0).$$ Clearly the difference between $\frac1x + \sin x$ and $\sin x$ will vanish in the limit $x\to\infty$, which implies that their roots will asymptotically agree, as claimed.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.