# How can I approximate an improper integral?

I have a function $f(x,y,z)$ such that
$\int_{R^3} f(x,y,z)dV$
is finite, and I want to approximate this integral.

I'm familiar with quadrature rules and monte carlo approximations of integrals, but I see some difficulties implementing them with in an infinite domain. In the monte carlo case, how does one go about sampling an infinite region (especially if the regions that contribute more significantly to the integral are unknown)? In the quadrature case, how do I find the optimal points? Should I simply fix an arbitrarily large region centered around the origin and apply sparse quadrature rules? How can I go about approximating this integral?

In one dimension, you can map your infinite interval to a finite interval using integration by substitution, e.g.

$$\int_a^b f(x)\,\mathrm dx \quad=\quad \int_{u^{-1}(a)}^{u^{-1}(b)}f(u(t))u'(t)\,\mathrm dt$$

Where $u(x)$ is some function that goes off to infinity in some finite range, e.g. $\tan(x)$:

$$\int_{-\infty}^{\infty}f(x)\,\mathrm dx \quad=\quad 2\int_{-\pi/2}^{\pi/2} f(\tan(t))\frac{1}{\cos(2t)+1} \,\mathrm dt$$

You can then use any regular numerical quadrature routine for the modified, finite integral.

Substitution for multiple variables is a bit trickier, but is quite well described here.

• That's very interesting... I never even considered the possibility of substitution! But does the choice of function $u(t)$ have any effect on the accuracy of the approximation? – Paul Jun 27 '12 at 13:50
• @Paul: Yes, definitely! The function $u(t)$ should be as smooth as possible such as to keep $f(u(t))$ as smooth as possible, thus allowing for a more accurate integration. – Pedro Jun 27 '12 at 14:15
• That's true, but what I had in mind was the rate at which u(t) converges to infinity? Does that also affect the accuracy? – Paul Jun 27 '12 at 14:36
• @Paul: I don't know if I understand your question correctly, but the function has to end up at infinity at one point or another. If it takes its time and then grows sharply, then this will introduce some large gradients in $f(u(t))$, which makes it more difficult to integrate and could thus affect accuracy. – Pedro Jun 27 '12 at 14:40
• Your derivative for the tangent was wrong; I fixed it. – J. M. Jun 30 '12 at 1:36

The standard way of doing it is to extract from the expression for $f(x)$ an exponential prefactor, transform that to $e^{-x^2}$, and then use Gaussian quadrature rules (or Gauss Kronrod) with this as a weight. If $f$ is smooth, this usually gives excellent results.

In $R^3$, the same works with weight $e^{-|x|^2}$, and appropriate cubature formulas can be found, e.g., in the book by Engels, numerical quadrature and cubature.

Online formulas are at http://nines.cs.kuleuven.be/ecf/

• That works well if your integrand is roughly exp(-x^2). If your integrand is approximately normal, but centered away from the origin, this approach can work poorly. – John D. Cook Jun 27 '12 at 15:58
• @JohnD.Cook: That's why I wrote ''extract an exponential prefactor, transform that to $e^{−x^2}$'', which usually involves a linear transformation, combining a translation moving the center to the origin, and rotations and scalings to make the level sets approximately spherical. The function itself can be quite far from normal. – Arnold Neumaier Jun 27 '12 at 16:28

For one-dimensional quadrature, you can check the book on Quadpack (a golden oldie but still very relevant in one-dimensional quadrature) and the techniques used in the algorithm QAGI, an automatic integrator for an infinite range.

Another technique is the double-exponential quadrature formula, nicely implemented for an infinite interval by Ooura.

For cubature, you can consult the Encyclopedia of cubature formulas by Ronald Cools.

• Note that double-exponential quadrature is a in essence a substitution method; you make a substitution that transforms your infinite-range integral into another infinite-range integral whose decay rate is, well, double-exponential... – J. M. Jun 30 '12 at 1:32
• @J.M. Correct. And you do it to get the best out of the Euler-Mclaurin summation formula for the trapezoid rule, as do the IMT transform and the TANH transform. A nice paper on the history of the DE written by one of the founding fathers can be found here – GertVdE Jun 30 '12 at 4:57

If you remember how quadrature worked, then you'll also know one way how to approximate infinite integrals. Namely: for quadrature, you approximate the function $f(x)$ you want to integrate by something similar, say a polynomial $\tilde f(x)$ (or a piecewise polynomial) for which you can write down the integral analytically. You get $\tilde f$ from $f$ by interpolation at interpolation points -- which will then be your quadrature points.

For infinite integrals, one approach uses the exact same thinking. For example, you could try to approximate $f(x)$ on the entire line using a simpler function, e.g. $\tilde f(x)=e^{-x^2}p(x)$ with a polynomial $p(x)$ that interpolates $f(x)e^{x^2}$ at a number of points. There are then simple formulas for computing the integral $\int_{-\infty}^{\infty} \tilde f(x) dx$. The choice of interpolation points follows similar logic as is done for the usual quadrature derivation.

If you want to use Monte Carlo integration, you could start by using importance sampling with a sampler that roughly approximates your integrand. The better your sampler matches your integrand, the less variance in your integral estimates. It doesn't matter than your domain is infinite as long as your sampler has the same domain.