pylab, integral of sinc function

[=?ISO-8859-1?Q?Sch=FCle_Daniel?=]

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.5484213527594 133e-08

ok, can be thought as zero

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

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.1404662440661 117

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?=]

my fault

In [31]: simple_integral (lambda x:sinc(x/pi), -1000, 1000)
Out[31]: 3.1404662440661 1

[Paul Rubin]

Schüle Daniel <uv**@rz.uni-karlsruhe.dewri tes:

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)) )

[=?ISO-8859-1?Q?Sch=FCle_Daniel?=]

[...]

>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 :)

[Paul Rubin]

Schüle Daniel <uv**@rz.uni-karlsruhe.dewri tes:

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.

[Robert Kern]

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.5484213527594 133e-08

ok, can be thought as zero

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

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.1404662440661 117

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.framewor k/Versions/2.5/lib/python2.5/site-packages/numpy-1.0.2.dev3521-py2.5-macosx-10.3-fat.egg/numpy/lib/function_base.p y
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

[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.5484213527594 133e-08

ok, can be thought as zero

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

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.1404662440661 117

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.com posite
S.integrate.Num pyTest S.integrate.cum trapz
S.integrate.__a ll__ S.integrate.dbl quad
S.integrate.__c lass__ S.integrate.fix ed_quad
S.integrate.__d elattr__ S.integrate.inf
S.integrate.__d ict__ S.integrate.new ton_cotes
S.integrate.__d oc__ S.integrate.ode
S.integrate.__f ile__ S.integrate.ode int
S.integrate.__g etattribute__ S.integrate.ode pack
S.integrate.__h ash__ S.integrate.qua d
S.integrate.__i nit__ S.integrate.qua d_explain
S.integrate.__n ame__ S.integrate.qua dpack
S.integrate.__n ew__ S.integrate.qua drature
S.integrate.__p ath__ S.integrate.rom b
S.integrate.__r educe__ S.integrate.rom berg
S.integrate.__r educe_ex__ S.integrate.sim ps
S.integrate.__r epr__ S.integrate.tes t
S.integrate.__s etattr__ S.integrate.tpl quad
S.integrate.__s tr__ S.integrate.tra pz
S.integrate._od epack S.integrate.vod e
S.integrate._qu adpack
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.qua d(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.124450933778 113, 6.8429604895257 158e-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.qua d(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.535454558897 3298, 1.4922039610659 907)

In [42]: S.integrate.qua d(lambda x:sinc(x/pi), -1000, 1000,limit=1000 )
Out[42]: (3.140466243937 5475, 4.5659823144674 379e-08)
Cheers,

f

ps - the 2nd number is the error estimate.