pylab, integral of sinc function

Discussion in 'Python' started by =?ISO-8859-1?Q?Sch=FCle_Daniel?=, Feb 19, 2007.

  1. Hello,

    In [19]: def simple_integral(func,a,b,dx = 0.001):
    ....: return sum(map(lambda x:dx*x, func(arange(a,b,dx))))
    ....:

    In [20]: simple_integral(sin, 0, 2*pi)
    Out[20]: -7.5484213527594133e-08

    ok, can be thought as zero

    In [21]: simple_integral(sinc, -1000, 1000)
    Out[21]: 0.99979735786416357

    hmm, it should be something around pi
    it is a way too far from it, even with a=-10000,b=10000

    In [22]: def ppp(x):
    ....: return sin(x)/x
    ....:

    In [23]: simple_integral(ppp, -1000, 1000)
    Out[23]: 3.1404662440661117

    nice

    is my sinc function in pylab broken?
    is there a better way to do numerical integration in pylab?

    Regards, Daniel
    =?ISO-8859-1?Q?Sch=FCle_Daniel?=, Feb 19, 2007
    #1
    1. Advertising

  2. my fault

    In [31]: simple_integral(lambda x:sinc(x/pi), -1000, 1000)
    Out[31]: 3.14046624406611
    =?ISO-8859-1?Q?Sch=FCle_Daniel?=, Feb 19, 2007
    #2
    1. Advertising

  3. [...]

    >> In [19]: def simple_integral(func,a,b,dx = 0.001):
    >> ....: return sum(map(lambda x:dx*x, func(arange(a,b,dx))))

    >
    > Do you mean
    >
    > def simple_integral(func,a,b,dx = 0.001):
    > return dx * sum(map(func, arange(a,b,dx)))
    >


    yes, this should be faster :)
    =?ISO-8859-1?Q?Sch=FCle_Daniel?=, Feb 19, 2007
    #3
  4. =?ISO-8859-1?Q?Sch=FCle_Daniel?=

    Paul Rubin Guest

    Schüle Daniel <-karlsruhe.de> writes:
    > In [19]: def simple_integral(func,a,b,dx = 0.001):
    > ....: return sum(map(lambda x:dx*x, func(arange(a,b,dx))))


    Do you mean

    def simple_integral(func,a,b,dx = 0.001):
    return dx * sum(map(func, arange(a,b,dx)))
    Paul Rubin, Feb 20, 2007
    #4
  5. =?ISO-8859-1?Q?Sch=FCle_Daniel?=

    Paul Rubin Guest

    Schüle Daniel <-karlsruhe.de> writes:
    > > return dx * sum(map(func, arange(a,b,dx)))

    > yes, this should be faster :)


    You should actually use itertools.imap instead of map, to avoid
    creating a big intermediate list. However I was mainly concerned that
    the original version might be incorrect. I don't use pylab and don't
    know what happens if you pass the output of arange to a function.
    I only guessed at what arange does.
    Paul Rubin, Feb 20, 2007
    #5
  6. =?ISO-8859-1?Q?Sch=FCle_Daniel?=

    Robert Kern Guest

    Schüle Daniel wrote:
    > Hello,
    >
    > In [19]: def simple_integral(func,a,b,dx = 0.001):
    > ....: return sum(map(lambda x:dx*x, func(arange(a,b,dx))))
    > ....:
    >
    > In [20]: simple_integral(sin, 0, 2*pi)
    > Out[20]: -7.5484213527594133e-08
    >
    > ok, can be thought as zero
    >
    > In [21]: simple_integral(sinc, -1000, 1000)
    > Out[21]: 0.99979735786416357
    >
    > hmm, it should be something around pi
    > it is a way too far from it, even with a=-10000,b=10000
    >
    > In [22]: def ppp(x):
    > ....: return sin(x)/x
    > ....:
    >
    > In [23]: simple_integral(ppp, -1000, 1000)
    > Out[23]: 3.1404662440661117
    >
    > nice
    >
    > is my sinc function in pylab broken?


    A couple things:

    1) The function is not from pylab, it is from numpy.

    2) Look at the docstring of the function, and you will notice that the
    convention that sinc() uses is different than what you think it is.

    In [3]: numpy.sinc?
    Type: function
    Base Class: <type 'function'>
    Namespace: Interactive
    File:
    /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/numpy-1.0.2.dev3521-py2.5-macosx-10.3-fat.egg/numpy/lib/function_base.py
    Definition: numpy.sinc(x)
    Docstring:
    sinc(x) returns sin(pi*x)/(pi*x) at all points of array x.

    --
    Robert Kern

    "I have come to believe that the whole world is an enigma, a harmless enigma
    that is made terrible by our own mad attempt to interpret it as though it had
    an underlying truth."
    -- Umberto Eco
    Robert Kern, Feb 20, 2007
    #6
  7. Schüle Daniel wrote:

    > Hello,
    >
    > In [19]: def simple_integral(func,a,b,dx = 0.001):
    > ....: return sum(map(lambda x:dx*x, func(arange(a,b,dx))))
    > ....:
    >
    > In [20]: simple_integral(sin, 0, 2*pi)
    > Out[20]: -7.5484213527594133e-08
    >
    > ok, can be thought as zero
    >
    > In [21]: simple_integral(sinc, -1000, 1000)
    > Out[21]: 0.99979735786416357
    >
    > hmm, it should be something around pi
    > it is a way too far from it, even with a=-10000,b=10000
    >
    > In [22]: def ppp(x):
    > ....: return sin(x)/x
    > ....:
    >
    > In [23]: simple_integral(ppp, -1000, 1000)
    > Out[23]: 3.1404662440661117
    >
    > nice
    >
    > is my sinc function in pylab broken?
    > is there a better way to do numerical integration in pylab?


    Pylab is mostly a plotting library, which happens (for historical reasons I
    won't go into) to expose a small set of numerical algorithms, most of them
    actually residing in Numpy. For a more extensive collection of scientific
    and numerical algorithms, you should look into using SciPy:

    In [34]: import scipy.integrate

    In [35]: import scipy as S

    In [36]: import scipy.integrate

    In [37]: S.integrate.
    S.integrate.Inf S.integrate.composite
    S.integrate.NumpyTest S.integrate.cumtrapz
    S.integrate.__all__ S.integrate.dblquad
    S.integrate.__class__ S.integrate.fixed_quad
    S.integrate.__delattr__ S.integrate.inf
    S.integrate.__dict__ S.integrate.newton_cotes
    S.integrate.__doc__ S.integrate.ode
    S.integrate.__file__ S.integrate.odeint
    S.integrate.__getattribute__ S.integrate.odepack
    S.integrate.__hash__ S.integrate.quad
    S.integrate.__init__ S.integrate.quad_explain
    S.integrate.__name__ S.integrate.quadpack
    S.integrate.__new__ S.integrate.quadrature
    S.integrate.__path__ S.integrate.romb
    S.integrate.__reduce__ S.integrate.romberg
    S.integrate.__reduce_ex__ S.integrate.simps
    S.integrate.__repr__ S.integrate.test
    S.integrate.__setattr__ S.integrate.tplquad
    S.integrate.__str__ S.integrate.trapz
    S.integrate._odepack S.integrate.vode
    S.integrate._quadpack


    These will provide dramatically faster performance, and far better
    algorithmic control, than the simple_integral:


    In [4]: time simple_integral(lambda x:sinc(x/pi), -100, 100)
    CPU times: user 7.08 s, sys: 0.42 s, total: 7.50 s
    Wall time: 7.58
    Out[4]: 3.1244509352

    In [40]: time S.integrate.quad(lambda x:sinc(x/pi), -100, 100)
    CPU times: user 0.05 s, sys: 0.00 s, total: 0.05 s
    Wall time: 0.06
    Out[40]: (3.124450933778113, 6.8429604895257158e-10)

    Note that I used only -100,100 as the limits so I didn't have to wait
    forever for simple_integral to finish.

    As you know, this is a nasty, highly oscillatory integral for which almost
    any 'black box' method will have problems, but at least scipy is nice
    enough to let you know:

    In [41]: S.integrate.quad(lambda x:sinc(x/pi), -1000, 1000)
    Warning: The maximum number of subdivisions (50) has been achieved.
    If increasing the limit yields no improvement it is advised to analyze
    the integrand in order to determine the difficulties. If the position of
    a
    local difficulty can be determined (singularity, discontinuity) one will
    probably gain from splitting up the interval and calling the integrator
    on the subranges. Perhaps a special-purpose integrator should be used.
    Out[41]: (3.5354545588973298, 1.4922039610659907)

    In [42]: S.integrate.quad(lambda x:sinc(x/pi), -1000, 1000,limit=1000)
    Out[42]: (3.1404662439375475, 4.5659823144674379e-08)


    Cheers,

    f

    ps - the 2nd number is the error estimate.
    Fernando Perez, Feb 23, 2007
    #7
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Replies:
    1
    Views:
    521
    Istvan Albert
    Jan 19, 2005
  2. Dr. Colombes
    Replies:
    3
    Views:
    658
    John Hunter
    Feb 23, 2005
  3. Replies:
    0
    Views:
    328
  4. Charles Krug
    Replies:
    2
    Views:
    5,341
    Charles Krug
    Apr 28, 2005
  5. =?ISO-8859-1?Q?Sch=FCle_Daniel?=

    pylab, matplotlib ... roots function question

    =?ISO-8859-1?Q?Sch=FCle_Daniel?=, Jan 21, 2007, in forum: Python
    Replies:
    2
    Views:
    337
    Robert Kern
    Jan 22, 2007
Loading...

Share This Page