[SciPy-user] how to make a sinc(x) function - "divide by zero" !

Sebastian Haase haase at msg.ucsf.edu
Mon Dec 12 12:30:08 CST 2005


Scott,
Thanks for your reply. (As Perry points out - my code actually works on 
numarray)
But in your code: Don't you get a "division by zero" in Numeric here ?
"where" is evaluation BOTH* results FIRST (over the ENTIRE array) and only 
AFTERWARDS "selects" based on the on the first argument
(* actually all "three" arrays)

Thanks,
Sebastian Haase


On Monday 12 December 2005 10:04, Scott Ransom wrote:
> Here is a simple version that I use.  It uses the Taylor
> expansion when the argument is near zero.  Note that this
> is using Numeric and umath right now.  This should be
> trivial to change to new scipy_core (which I am not using yet):
>
> def sinc(xs):
>     """
>     sinc(xs):
>         Return the sinc function [i.e. sin(pi * xs)/(pi * xs)]
>             for the values xs.
>     """
>     pxs = umath.pi*xs
>     return Numeric.where(umath.fabs(pxs)<1e-3,
>                          1.0-pxs*pxs/6.0, umath.sin(pxs)/pxs)
>
> Scott
>
> On Monday 12 December 2005 12:39 pm, Arnd Baecker wrote:
> > On Mon, 12 Dec 2005, Sebastian Haase wrote:
> > > Hi,
> > > I was trying to implement a "sinc" [sin(x)/x] in numarray. But the
> > > "where"-semantics makes it choke on the x=0: 0/0 case.
> > > (This should of course work for x being an array - so "if" is no
> > > option) The best I could come up with was:
> > >
> > > def sinc(r):
> > >     na.Error.pushMode(all="ignore")
> > >     a = na.where(r, na.divide(na.sin(r),r), 1)
> > >     na.Error.popMode()
> > >     return a
> > >
> > > but I still seem to get  a warning...
> > >
> > > >>> F.sinc(0)
> > >
> > > Warning: Encountered invalid numeric result(s)  in divide
> > > 1.0
> > >
> > > What is a better way of doing this ?
> >
> > Yuo could try to use `vectorize` (warning: untested code)
> >
> > def sinc(x):
> >     if x==0.0:                # presumably better to check for small
> > x return 0.0            # here ...
> >     else:
> >        return sin(x)/x
> >
> > sinc_vectorized=scipy.vectorize(sinc)
> >
> > HTH,
> >
> > Arnd
> >
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user at scipy.net
> > http://www.scipy.net/mailman/listinfo/scipy-user



More information about the SciPy-user mailing list