[SciPy-User] integrals and quad
nicky van foreest
Thu Aug 16 03:13:42 CDT 2012
Hi Josef and Sasha,
Thanks for the hints. I tried cumtrapz, and it gives nice results, at
least for the trivial integration of f(x) = x.
Now that I know cumtrapz, I realize I could have found this myself,
but just reading the scipy.integrate docs...
On 15 August 2012 22:36, Oleksandr Huziy <email@example.com> wrote:
> I am not sure if it falls under the elegant, but what if you take the
> integrals for parts and then apply cumsum:
>>>> import numpy as np
>>>> t1 = np.linspace(0,3,50)
>>>> t2 = t1[1:]
>>>> t1 = t1[:-1]
>>>> parts = map(lambda x1,x2: quad(f,x1,x2), t1,t2)
>>>> F = np.cumsum(parts)
> Probably you'd want to add 0 at the start of the result...
> Oleksandr (Sasha) Huziy
> 2012/8/15 <firstname.lastname@example.org>
>> On Wed, Aug 15, 2012 at 4:15 PM, nicky van foreest <email@example.com>
>> > 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)) # must be wasteful
>> > t = linspace(0,3, 50)
>> > FF = cumsum(f(t))*(t-t) # simple, but inaccurate, note that
>> > t-t 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?
>> cumtrapz is the only one that works (when I looked at this)
>> I also tried odeint for this once before, but didn't really use it.
>> > thanks
>> > Nicky
>> > _______________________________________________
>> > SciPy-User mailing list
>> > SciPy-User@scipy.org
>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>> SciPy-User mailing list
> SciPy-User mailing list
More information about the SciPy-User