[SciPy-user] Integrating Array

Anne Archibald peridot.faceted@gmail....
Thu Oct 18 12:15:26 CDT 2007


On 18/10/2007, rgold@lsw.uni-heidelberg.de <rgold@lsw.uni-heidelberg.de> wrote:

> I need to integrate a set of data points stored in a 1-dim array. That
> means I can only use scipy.integrate.simps and scipy.integrate.trapz but
> these methods are not as accurate as the routines for function-input like
> scipy.integrate.quadrature which efficiently use more samples where the
> integrand oscillates faster,etc.
> The problem is that I need a high accuracy because the integrand IS
> oscillating fast!
>
> I already interpolated the data using the scipy.interpolate.splrep but
> interestingly it turns out that simpsons rule now produces NANs and I am
> left with good old trapz!
>
> Is there a way of implementing quadrature for data-like-arrays?
> Still I could try romberg. But then I have to interpolate exactly in such
> a way that the "number of samples = positive power of 2) +1" (equally
> spaced if I remember correctly).
>
> Is this the best way or does anyone have another idea?

The reason the clever routines (e.g., scipy.integrate.quad) are able
to produce high accuracy is because they are able to sample your
function more densely in regions where it is more complicated. Between
samples, they generally make some very simple assumption about the
function's behaviour - a polynomial of low order, usually. If you have
already sampled your function and are not willing to sample it
further, there's not much that can be done to improve on simps. If
your function were very smooth, you could try very high order
polynomials, but this  is unlikely to help you much. The real problem
is that you only know your function at some points, and interpolation
is the best guess you have at what happens between them. The short
answer is if your function oscillates between points, you need more
points.

The best and easiest choice is to simply use quad on your function,
and let it decide where to sample. From your question, you presumably
have a good reason not to do that. If it's that your evaluation
routine is much more efficient when you give it a whole vector to
evaluate at once, there are choices:

* scipy.integrate.quadrature is based on Gaussian quadrature. The idea
here is to evaluate your function at carefully-chosen abscissas, then
use carefully-chosen weights to produce a quadrature rule. The
abscissas and weights can be chosen so that you get exact integrals
for polynomials up to a certain order times a weight function. Anyway,
it's good, and it evaluates your function at a vector of points at
once. scipy.integrate.quadrature evaluates it on denser and denser
grids until the result appears to converge. (If you need special
calling conventions for your function, you can look into the code of
scipy.integrate.quadrature, it's not very complicated.)

* If your function can only be evaluated at an evenly-spaced grid, you
can use scipy.integrate.simps on finer and finer grids until it
appears to have converged. If you have slightly more flexibility, you
can subdivide individual grid steps separately (so that regions where
the function is more complicated get evaluated more frequently).

* On the remote chance that you are getting Fourier coefficients and
wanting to integrate the time-domain function, some clever algebra
involving the inverse Fourier transform lets you produce integrals
directly (using an IFFT if desired).

We may be more able to help you if you tell us more about your problem
(for example why you're not using quad).

Good luck,
Anne


More information about the SciPy-user mailing list