Dear All
Here is the long overdue summary of the responses to my query concerning
numerical integration. Thank you so much to all those who responded.
Kind regards
Lars Thomassen
Presumably, if you require the integral to a certain degree of accuracy, you
would be able to obtain estimates by using intervals of size [0,10^n] for
n=1,2,3,.... and stopping when you reach the required level of accuracy. If
the integral exists on [0,+infinity) then you should have the required
convergence, so whilst it is not the most analytically sound way of
producing the integral it should give a suitable approximation.
-------
See Bernsten, Espelid and Genz: Algorithm 698: Dcuhre: An adaptive
multidemensional integration routine for a vector of integrals. Transactions of
Mathematical Software, 17(4):452-456, 1991.
This algorithm allows for limits of +/- infinity. I downloaded a version in C
from Genz's website at Washington State University a few years ago. He has also
done other papers on the subject.
----
This is a guess.
To integrate f() find two analytically integrable functions g() and h()
such that g <= f <= h on the interval [a, +infinity).
The tail area of f() is bounded below by the tail area of g() and bounded
above by the tail area of h().
Appropriate choices of g, h, and the tail [a, +infinity) allow you to
choose the accuracy of the integral.
----
Try using Gauss Laguerre rule.
See
Abramowitz and Stegun
Handbook of mathematical functions.
Dover
ISBN 486-61272-4
-----
Make sure the function you want to integrate is well-behaved near
infinity. Can it be approximated by for instance a negative exponential
or a power of x less than -1. If so, the integral is finite and from a
certain point outwards the contribution to the numerical integration can
be calculated. If it is not, the integral diverges.
-----
One of the best ways is to use a 1-1
mapping so that if your variable of integration x is defined on [0,inf) then
(tan inverse x) will be defined on [0, Pi/2).
-----
One method is to split the integral into one from 0 to 1 and one from 1 to
infinity, then transform the latter integral so that you integrate with
respect to 1/x, so that the infinity limit becomes 1.
----
Romberg integration will do it. I have it implemented (vectorized) as
a function in R available at
www.luc.ac.be/~jlindsey/rcode.html
------
Two very good references are: Thisted (1988) Elements of statistical
computing and the "Numerical Recipes in C". The fundamental trick is to
apply a transformation to the integrand, such that the interval becomes
finite (for example exp(-x)). Thisted has a nice discussion on how to
choose such a transformation.
------
Look up 'Gaussian quadrature' in a standard numerical analysis text,
or look up Splus on integration.
------
I'd recommend "Numerical Methods that Work" by Forman S. Acton. His
view is that one should be prepared to use several different ideas
(which he describes) and choose the most appropriate one for the
integral in question.
------
Ronald Thisted, "Elements of Statistical Computing"
(Chapman and Hall, 1988), on numerical integration. It gives a
number of techniques suited to different integrands, and gives some
references as well.
----
NAG have a comprehensive library of subroutines, including a large selection
of subroutines on numerical integration (quadrature). Two of these deal
with infinite, or semi-infinite, ranges. You can get a full description of
the two subroutines, complete with references on the Web. The Web address
below will get you to the appropriate page:-
http://www.nag.co.uk/numeric/FL/manual/html/genint/FLlibconts.asp
-----
|