# specific improper integral involving erf

I have encountered an integral, and kindly ask for help with a solution. It is beyond my own capabilities, and neither Maple nor Mathematica were of any help: $$\int_{1}^{\infty} \left[\mathrm{erf}\left(px\right)\right]^{\alpha} \exp\left(-p^2x^2\right)\,\frac{\mathrm{d}p}{p},$$ where $\alpha$ is an integer $> 0$, and $x>0$.

If the lower integration boundary were $0$ instead of $1$, a Laplace transform might be of help (with the transformation $p\rightarrow k=p^2x^2$). However, integrating that part separately leaves me with another integral I can't tackle (with integration boundaries $0$ and $x^2$, or $0$ and $1$ before the transformation).

I would also be grateful for hints for efficiently solving this integral numerically if it cannot be simplified further. Again, for $0$ as the lower integration boundary, (generalized) Gauss-Laguerre quadrature would help for one part, but the other is still a problem (presumably Gauss-Legendre quadrature is not conducive because of the singularity at $0$).

• Mathematica's NIntegrate works just fine for a numerical evaluation (no hope for a closed form solution, of course). – Carlo Beenakker Sep 24 '13 at 12:22
• NIntegrate does work, of course. However, it is not very useful since this integral is part of the argument of a Fourier sum (over x) with up to several million terms. I am currently solving it "brute force" with the midpoint rule in a C code (on a logarithmic grid), but was hoping that there is an ingenious way to maybe find a special function somewhere in this integral (if there's no closed-form solution), or another clever way to make evaluation more efficient. – Monolithus Sep 24 '13 at 12:47

Here's some ideas. Your integrand is non-oscillatory with good decay so you should only need <100 quadrature pts for machine accuracy if you use a high-order method. I suggest plain trapezoid rule with Alpert end correction at the bottom end and truncation in the Gaussian tail at top end. Eg translate code from http://mpspack.googlecode.com/svn/trunk/@quadr/QuadNodesInterval.m http://mpspack.googlecode.com/svn/trunk/@quadr/QuadSmoothExtraPtNodes.m If erf is slow build a look-up table using piecewise polynomial interpolation on Chebyshev grids.

Also: sub y=px and your function of $x$ is an antiderivative of $(\mbox{erf} y)^\alpha e^{-y^2} / y$; if this is to be summed over $x$ as you say, compute the antiderivative together as you go along the $x$ integral.

Also: try to find recurrence relation for $\alpha$ to $\alpha+1$ if you care about many alpha cases (I didn't spot anything yet).

• thanks for your ideas. Only one (fixed) $\alpha$ is needed, most likely $\alpha=1$ or $\alpha=2$. – Monolithus Sep 25 '13 at 6:42

As others have pointed out, there's no hope for a closed-form solution. To evaluate it numerically in an efficient way, I would first rewrite this as $$\int_x^\infty (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y}\;.$$ Now if $y$ is larger than about $10$ or so, then the integrand is equal to $0$ for all practical purposes, so that the integral is pretty well approximated by $$\int_x^{10} (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y}\;,$$ which can easily be evaluated to machine precision by any quadrature formula (provided you can evaluate the error function). If you want your code to be fast, just precompute the integrand at a few thousand grid points...

If you are interested in small relative error rather than just small absolute error, then you can also get asymptotics for large $x$. To leading order, you'll get a pretty good approximation by $$\int_x^\infty (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y} \approx {1\over x} \int_x^\infty \exp(-y^2) dy = {\sqrt{\pi} \, \mathrm{erfc}(x) \over 2x} \approx {e^{-x^2} \over 2x^2}\;.$$

• Thanks for your suggestions. The error function is easily computed, for instance with the boost library boost, albeit for C++ (which is no problem). How did you do the first step in the large-x asymptotic approximation? – Monolithus Sep 25 '13 at 6:49
• Replace $\mathrm{erf}(y)$ by $1$ because $y$ is large and replace $1/y$ by $1/x$ because $\exp(-y^2)$ is going to be concentrated near $x$. – Martin Hairer Sep 27 '13 at 14:18