pylab, integral of sinc function

[Fernando Perez]

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.