Re: pylab, integral of sinc function
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=-1,b=1 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. -- http://mail.python.org/mailman/listinfo/python-list
pylab, integral of sinc function
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=-1,b=1 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 -- http://mail.python.org/mailman/listinfo/python-list
Re: pylab, integral of sinc function
my fault In [31]: simple_integral(lambda x:sinc(x/pi), -1000, 1000) Out[31]: 3.14046624406611 -- http://mail.python.org/mailman/listinfo/python-list
Re: pylab, integral of sinc function
Schüle Daniel [EMAIL PROTECTED] 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))) -- http://mail.python.org/mailman/listinfo/python-list
Re: pylab, integral of sinc function
[...] 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 :) -- http://mail.python.org/mailman/listinfo/python-list
Re: pylab, integral of sinc function
Schüle Daniel [EMAIL PROTECTED] 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. -- http://mail.python.org/mailman/listinfo/python-list
Re: pylab, integral of sinc function
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=-1,b=1 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 -- http://mail.python.org/mailman/listinfo/python-list
