1
$\begingroup$

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$).

Improve this question
$\endgroup$
2
  • $\begingroup$ Mathematica's NIntegrate works just fine for a numerical evaluation (no hope for a closed form solution, of course). $\endgroup$
    – Carlo Beenakker
    Commented Sep 24, 2013 at 12:22
  • $\begingroup$ 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. $\endgroup$
    – Monolithus
    Commented Sep 24, 2013 at 12:47

2 Answers 2

Reset to default
2
$\begingroup$

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).

Improve this answer
$\endgroup$
1
  • $\begingroup$ thanks for your ideas. Only one (fixed) $\alpha$ is needed, most likely $\alpha=1$ or $\alpha=2$. $\endgroup$
    – Monolithus
    Commented Sep 25, 2013 at 6:42
1
$\begingroup$

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}\;. $$

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Monolithus
    Commented Sep 25, 2013 at 6:49
  • $\begingroup$ 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$. $\endgroup$
    – Martin Hairer
    Commented Sep 27, 2013 at 14:18

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .