Deal with singularities for numerical integration 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)?
 A: 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.
A: 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.
A: $\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));
}

