[SciPy-User] integrals and quad

[nicky van foreest]

Hi,

Given some function f it is easy with scipy.integrate.quad to compute
the integral of f for some given endpoint. However, I need the
integral at many endpoints, that is, I want to plot \int_0^t f(x) dx.
How can this be done in an efficient and elegant way?

To illustrate I used the following code.

from numpy import cumsum, linspace, vectorize
from scipy.integrate import quad
from  pylab import plot, show

def f(x):
    return x

F = vectorize(lambda t: quad(f, 0, t)[0]) # must be wasteful


t = linspace(0,3, 50)

FF = cumsum(f(t))*(t[1]-t[0])  # simple, but inaccurate, note that
t[1]-t[0] is the grid size, a bit like dx in the integral

plot(t,f(t))
plot(t, F(t))
plot(t, F)
show()

I suspect that calling F at many values is wasteful, since the
integral is evaluated at the same points many times. The trick with
using cumsum must save some work (an O(n) algo), but is less accurate
as is shown by the graphs. So, I don't like to use cumsum, and I also
don't like to use a vectorized quad. Is there something better?

thanks

Nicky