Algorithm for calculating definite integrals with bounds at infinity

Question

Suppose I have an integral that's bounded on one (or both) end by (-)infinity. AFAICT, I can't analytically solve this problem, it takes brute force (e.g. using a Left Riemann Sum). I'm having trouble generalizing the algorithm so that it sets the proper subdivisions; I'll either do far too much work to calculate something trivial, or not do nearly enough and have huge aliasing errors.

Answering in any language is cool, but maybe someone with better google-fu can end this quickly. :)

Is what I'm looking for as impossible as trying to measure the British coastline?

Answer 1

If you're looking for an exhaustive solution, you'll probably need a symbolic math package of some kind. Reasons:

- If you're indeed looking for solutions to indefinite integrals (ones with no bounds at all, not just finite ones) then I don't know any other way except doing the digital equivalent of a lookup table or a symbolic algorithm (like knowing the integral of x^n dx is x^(n+1) / (n+1), plus an arbitrary expression that doesn't depend on x).

• You can get somewhere with numeric integration on definite integrals. If one or both of the limits is +/- infinity, then you can probably approximate this with finite numbers that are "infinite enough," but then again, that gets you only so far as well, because what about things like the integral of 1/x from 1 to infinity? The answer is infinity, but it diverges as ln(x). There is no value that's "infinite enough" to give you the answer. Which leads back to a symbolic solution.

Answer 2

Use a change of variable. For example, x -> 1/x and note that the integral from a to b of f(x) dx is equal to the integral from 1/b to 1/a of (1/x^2)f(1/x) dx. Another handy one is the change of variable x -> -log(x), where the integral from a to infinity of f(x) dx is equal to the integral from 0 to e^(-a) of f(-log(x))/x dx.

Various open source packages contain routines to do this sort of thing. Someone else will have to recommend one, because I wrote my own when open source for math functions was non-existent.

The Numerical Recipes in C book is available online http://www.nrbook.com/a/bookcpdf.php You may be able to cobble something together from that. The chapter you are interested in is number 4. Be aware that array indexes are based off 1 rather than 0. Some of the algorithms could stand improvement, but I'll not go into that here.

Answer 3

There are several ways to proceed, most of which involve trying to understand the behavior of your integrand. Often, there is a transformation x -> z(x) with a finite z(infinity) so you can transform your unbounded integral to a bounded one.

Also, it is often the case that you can analyze the "asymptotic" behaviour of the integrand as x goes to + and - infinity, so you can at least approximately work out the contributions from x>x+ and x< x- and then just do the definite integral between x- and x+.

There are plenty of good books on numerical integration. One that's used a lot in the physical sciences is Numerical Recipes. (Although you should rarely use the code directly!)