Estimating the value of an improper integral numerically

Question

My question is how can I estimate the value of an improper integral from $[0,\infty)$ if I only have a programming routine that gives me the function evaluated at 100 data points, or 100 values of $x$. I don't actually have the expression for the function at all. Assuming that the function is continuous on $[0,\infty)$, how do I go about doing this problem. do I possibly use Riemann sums? Note that my $f(x)$ is actually a matrix function, meaning that $f(x)$ actually is a matrix with elements that are functions of $x$. Also, this will be implemented in Python (a programming language) so computational methods would be the best.

NOTE: for those who know Python, I have solved a matrix differential equation using scipy.integrate.odeint which gives me numerical results of $f(x)$ evaluated at 100 data points. The routine I am talking about is one that runs this function for different inputs for x. EDIT:

I need to add two things:

1. As I mentioned, my integrand is continuous from 0 to infinity. I just don't know the functional for, of my integrand. I can get the values of $f(x)$ for any x I choose

2. I don't want to interpolate or do regression

Answer 1

You can transform the integral by substituting some function for $x$ that goes to infinity at a finite value, e.g. $\tan u$. That yields

$$ \int_0^\infty f(x)\,\mathrm dx=\int_0^{\pi/2}\frac{f(\tan u)}{\cos^2 u}\,\mathrm du\;. $$

This integral you can evaluate using standard techniques. Of course this only works if your integral converges, or equivalently, if $f(\tan u)\rightarrow0$ sufficiently quickly as $\cos^2u\rightarrow0$ at $u=\pi\,/\,2$.

Edit: We clarified in the chat that the problem had been incompletely stated. The function can indeed be evaluated at any selected $x$ value, but this is very expensive and automatically yields function evaluations at $100$ evenly spaced values of $x$. In this case, standard (e.g. Gaussian) quadrature in $u$ makes no sense, since it's better to use $100$ points for the trapezoidal rule than $1$ point for Gaussian quadrature.

This of course kills much of the value of transforming to $u$, but two potential advantages remain. First, the linear approximation implied in the trapezoidal rule may be better in the transformed function (or it may not, depending on $f$). Second, you can handle the contribution at $x=\infty$ gracefully, whereas in the improper integral there's no trapezoid you could construct with one side at $x=\infty$.

This is what the computation would look like concretely (with $x_i$ the evenly spaced values of $x$):

\begin{eqnarray} &&\int_0^\infty f(x)\,\mathrm dx=\int_0^{\pi/2}\frac{f(\tan u)}{\cos^2 u}\,\mathrm du\\ &&\simeq\frac12\left(\arctan x_1f(0)+(\pi/2-\arctan x_{100})\lim_{u\to\infty}\frac{f(\tan u)}{\cos^2 u}+\sum_{i=1}^{100}\frac{\left(\arctan x_{i+1}-\arctan x_{i-1}\right)f\left(x_i\right)}{\cos^2\arctan x_i}\right)\\ &=&\frac12\left(\arctan x_1f(0)+(\pi/2-\arctan x_{100})\lim_{u\to\infty}\frac{f(\tan u)}{\cos^2 u}+\sum_{i=1}^{100}\left(\arctan x_{i+1}-\arctan x_{i-1}\right)f\left(x_i\right)\left(1+x_i^2\right)\right)\;. \end{eqnarray}

This requires you to know the limit of $f(\tan u)\,/\cos^2 u$. If you can't determine that, the next best thing is to use a rectangle instead of a trapezoid at that end, but that would make the error closer to the one incurred by Steven's method, so in that case you might be better off using that method directly. The best way to tell is to try it out and see how rapidly the results converge in each case when you increase the number of points.

Answer 2

You are going to have to assume that your last point, $b$, is such that $\int_{b}^{\infty}f(x)\,dx$ is neglegible. If that is really true, then a Reimann sum approximation of $\int_a^b f(x)\,dx$ will do fine. Since you seem to have little control over the value of $\Delta x$, I wouldn't bother getting more creative than that.

You should at least plot your data to check that the transient assumption is valid. If it's not, then you need to figure out some strategy to make up for the missing area.

Computing $\int_a^b f(x)\,dx$ depends a lot on how much you know about $f(x)$. You seem to be able to calculate $f(x)$ but, if don't want to interpolate, you may as well treat it as raw data.

The book Numerical Recipes In C can be found online here and probably many other places. Don't let the C part bother you. It is very readable and a good source of ideas and methods.

I have two copies of Numerical Methods That Work by Forman S. Acton. The book is badly bound and tends to fall apart. It is very practical but a bit more work to read. It has made me look smarter than I am many times.