Deal with singularities for numerical integration

Question

I am trying to find a good way to numerically integrate $$\int_0^\infty \frac{e^{ -x^2/(2 a^2)} }{a \sqrt{2 \pi}} \frac{2 \log( a^2 ) }{\pi^2 (a^2-1)} da$$ for $x \in (-\infty, \infty)$. I believe there is no closed form analytical solution for this integral. In the end I want to plot this integral as a function of x.

My main issue with using some simple (adaptive) quadrature procedures is the singularities at 0 (due to the $\log a^2$ term) and at 1 (due to the $1/(a^2-1)$ term) - plus the integral being to $\infty$. I would assume that one could find some good substitution that might remove both singularities and turn this into an integral on finite intervals, but I lack the experience or intuition about what a good substitution would be in this case.

I thought of $u=\cos(a)$ or $u = \log a^2$, but these do not do the job. Or perhaps I could split the integral into e.g. 0 to 0.5, 0.5 to 2 and 2 to $\infty$. For the first of these $u=\log a^2$ does remove the singularity. Then I need some transformation that removes the signularity in the second integral and one that turns the final interval into a finite one. Are there obvious ones I should consider.

Or is there some more obvious or simple way to do this numerical integration (particularly, with a method that is implemented in some R package)?

Answer 1

Let $F(x)=\int_0^\infty \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da$. We remark that $F(0)$ fails to exist.

Next, we can write for $x\ne 0$

$$\begin{align} F(x)&=\int_0^\infty \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da\\\\ &=\int_0^1 \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da+\int_1^\infty \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da\\\\ &=\int_0^1 \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da+\int_0^1 \frac{e^{-x^2a^2/2}}{1/a}\frac{\log(1/a^2)}{1/a^2-1}\frac{1}{a^2}\,da\\\\ &=\int_0^1 \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da+\int_0^1 ae^{-x^2a^2/2}\frac{\log(a^2)}{a^2-1}\,da \tag 1 \end{align}$$

Note that $\frac{\log(a^2)}{a^2-1}=\frac{2}{a+1}\frac{\log(a)}{a-1}$ is well-behaved near $a=1$ since $\log(a)=(a-1)+O(a-1)^2$ for $a$ near $1$. Hence, both integrals on the right-hand side of $(1)$ are well-behaved near $a=1$

The second integral on the right-hand side of $(1)$ is also well-behaved near $a=0$ since $\lim_{a\to 0}a\log(a)=0$.

Finally, we see that for fixed $x\ne 0$, $\lim_{a\to 0} \frac{e^{-x^2/(2a^2)}}{a}\,\log(a^2)=0$ and hence the integral on the right-hand side of $(1)$ is also well-behaved near $a=0$.

Hence, for numerical evaluation, we need only redefine the integrand as to account for the removable discontinuities.

---

All of that said, there is still a potentially challenging issue left to resolve. For $x=0$, the first integral on the right-hand side of $(1)$ diverges. We can handle this issue as follows.

To facilitate numerical analysis, we proceed by integrating by parts with $u=\frac{a^2\log(a^2)}{a^2-1}$ and $dv=\frac{e^{-x^2/(2a^2)}}{a^3}$. Then, we find that

$$\begin{align} \int_0^1 \frac{e^{-x^2/(2a^2)}}{a}\,\frac{\log(a^2)}{a^2-1}\,da&=\frac{e^{-x^2/2}}{x^2}-\frac{1}{x^2}\int_0^1 e^{-x^2/(2a^2)}\left(\frac{2a\log(a^2)}{a^2-1}-\frac{2a^3\log(a^2)}{(a^2-1)^2}+\frac{2a}{a^2-1} \right)\,da \tag 2 \end{align}$$

The integral on the right-hand side of $(2)$ is well-behaved for $x=0$. In fact, at $x=0$ the integral is equal to $1$.

Finally, we write

$$\begin{align} F(x)&=\frac{e^{-x^2/2}}{x^2}-\frac{1}{x^2}\int_0^1 e^{-x^2/(2a^2)}\left(\frac{2a\log(a^2)}{a^2-1}-\frac{2a^3\log(a^2)}{(a^2-1)^2}+\frac{2a}{a^2-1} \right)\,da\\\\ &+\int_0^1 ae^{-x^2a^2/2}\frac{\log(a^2)}{a^2-1}\,da \tag 3 \end{align}$$

where all of the integrals in $(3)$ have removable discontinuities.

Answer 2

I would tackle this by substituting $a = \frac{x}{\sqrt{2y}}$ making the integral an exponential of fixed width multiplied by a function that is algebraic up to a logarithmic term, the $x$ dependence is then in this latter factor. This then makes Gauss–Laguerre quadrature quite effective. The substitution leads to the integral:

$$\frac{1}{\sqrt{2\pi}\pi^2}\int_0^{\infty}\frac{\exp(-y)}{y-\frac{x^2}{2}}\log\left(\frac{2y}{x^2}\right)dy$$

Then we note that the integrand has an integrable singularity at $y = 0$ and has a removable singularity $y = \frac{x^2}{2}$, see Dr. MV's solution for details. Due to the logarithmic term, one can set up a generalized Gauss-Laguerre quadrature with weight function an exponential multiplied by a logarithm, see here for details.

Answer 3

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ You can redefine the term $\ds{{\ln\pars{a^{2}} \over a^{2} - 1}}$. Lets $\mc{MP}$ the $machine\ precision$ \begin{align} {\ln\pars{a^{2}} \over a^{2} - 1} & = {\ln\pars{1 + \bracks{a^{2} - 1}} \over a^{2} - 1} \approx {\pars{a^{2} - 1} - \pars{a^{2} - 1}^{2}/2 + \pars{a^{2} - 1}^{3}/3 - \pars{a^{2} - 1}^{4}/4 \over a^{2} - 1} \\[5mm] & = 1 - {1 \over 2}\pars{a^{2} - 1} + {1 \over 3}\pars{a^{2} - 1}^{2} - {1 \over 4}\,\pars{a^{2} - 1}^{3} \end{align} When $\ds{\verts{{1 \over 4}\,\pars{a^{2} - 1}^{3}} < \mc{MP} \implies \verts{a^{2} - 1} < \pars{4\mc{MP}}^{1/3}}$, you can use the $approximation$ \begin{align} {\ln\pars{a^{2}} \over a^{2} - 1} & \approx 1 - {1 \over 2}\pars{a^{2} - 1} + {1 \over 3}\pars{a^{2} - 1}^{2} = 1 - {1 \over 2}\pars{a^{2} - 1}\bracks{1 - {2 \over 3}\pars{a^{2} - 1}} \end{align}

---

The following $\texttt{C++}$ code implements the above argument:

#include <cmath>
#include <cfloat> // DBL_EPSILON is the Machine Precision for the type double 
const double FACTOR23 = 2.0/3.0,TOL = 1.05*pow(4.0*DBL_EPSILON,1.0/3.0); 
// TOL includes a 5 % additional factor to avoid rounding errors.

double f(double a)
{
 a = a*a - 1.0;
 return (abs(a) > TOL) ? (log(1.0 + a)/a):(1.0 - 0.5*a*(1.0 - FACTOR23*a));
}