Opened 12 years ago

Closed 12 years ago

#10239 closed defect (fixed)

Bessel method forces bessel_Y to use pari algorithm

Reported by: jvkersch Owned by: jason, jkantor
Priority: major Milestone: sage-4.6.1
Component: numerical Keywords: functions, bessel, pari
Cc: Merged in: sage-4.6.1.alpha2
Authors: Jeroen Demeyer Reviewers: Joris Vankerschaver
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description

When computing the Bessel Y-function, the Bessel function calls bessel_Y with a default argument of algorithm='pari'. However, bessel_Y throws an exception when using anything else than scipy or maxima.

sage: f = Bessel(0, typ='Y')
sage: f(1.5)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)

/Users/Joris/Desktop/sage/build/sage-4.6.alpha2/<ipython console> in <module>()

/Users/Joris/Desktop/sage/build/sage-4.6.alpha2/local/lib/python2.6/site-packages/sage/functions/special.pyc in __call__(self, z)
   1084             return bessel_K(nu,z,algorithm=s,prec=p)
   1085         if t == "Y":
-> 1086             return bessel_Y(nu,z,algorithm=s,prec=p)
   1087         
   1088     def plot(self,a,b):

/Users/Joris/Desktop/sage/build/sage-4.6.alpha2/local/lib/python2.6/site-packages/sage/functions/special.pyc in bessel_Y(nu, z, algorithm, prec)
    938         return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    939     else:
--> 940         raise ValueError, "unknown algorithm '%s'"%algorithm
    941     
    942 class Bessel():

ValueError: unknown algorithm 'pari'

The following works fine:

sage: f = Bessel(0, typ='Y', algorithm='scipy')
sage: f(1.5)
0.382448923798000
sage: g = Bessel(0, typ='Y', algorithm='maxima')
sage: g(1.5)
0.382448923797759

As a first solution, it would be easy to change the calling convention of bessel_Y to ensure that maxima or scipy is used. However, it seems like pari does a better job at computing the Bessel functions, and therefore it might be better in the long run to extend bessel_Y to use pari (if at all possible).

Attachments (2)

10239_bessel_Y.patch (8.8 KB) - added by jdemeyer 12 years ago.
10239_bessel_Y_reviewer.patch (8.7 KB) - added by jvkersch 12 years ago.

Download all attachments as: .zip

Change History (10)

comment:1 Changed 12 years ago by jdemeyer

  • Component changed from misc to numerical
  • Milestone set to sage-4.6.1
  • Owner changed from jason to jason, jkantor

The problem is that PARI has no function to compute the Y Bessel function. It seems that PARI can compute the following Bessel functions:

besselh1  besselh2  besseli   besselj   besseljh  besselk   besseln

I can make a patch which simply changes the default algorithm to "maxima" for bessel_Y.

Changed 12 years ago by jdemeyer

comment:2 Changed 12 years ago by jvkersch

Too bad that there's no pari implementation: for Bessel-J it is just so fast:

sage: timeit(r"bessel_J(1, 2.3, algorithm='pari')")
625 loops, best of 3: 110 µs per loop
sage: timeit(r"bessel_J(1, 2.3, algorithm='scipy')")
625 loops, best of 3: 480 µs per loop
sage: timeit(r"bessel_J(1, 2.3, algorithm='maxima')")
125 loops, best of 3: 5.27 ms per loop

Still, using the maxima implementation would be good for now, as the docs mention that there are some accuracy issues with scipy. Also, thanks for your quick response -- it feels like I never left Ghent and we are still talking about Sage in the backyard of the department!

comment:3 Changed 12 years ago by jvkersch

  • Status changed from new to needs_work

The patch seems to work well. The Bessel method now behaves as intended, and the patch also addresses some other inconsistencies in the code (checking the precision for the scipy code, and some documentation issues).

There is one doctest that fails on my machine (64bit macbook pro, sage 4.6):

sage: Bessel(6, algorithm="scipy")(pi)
Expected:
    0.0145459669825000 + 5.34392507678000e-18*I
Got:
    0.0145459669825000 + 5.34410157169000e-18*I

This looks like inevitable round-off error in the complex part. Can we hit this doctest with the #random tag or is there a more subtle way of dealing with this? IMHO, this is also a bug in scipy: for real input, the Bessel functions of the first and second kind should be exactly real, not up to some small complex drift. This will cause all sorts of problems with plotting, etc.

comment:4 Changed 12 years ago by jvkersch

  • Status changed from needs_work to needs_review

The problem with scipy is the way in which the argument is presented to the Bessel function: when the argument is real, scipy returns a real value for the Bessel function:

sage: import scipy.special
sage: scipy.special.jv(6, complex(float(pi)))
(0.014545966982505558+5.3441015716885544e-18j)
sage: scipy.special.jv(6, float(pi))
0.014545966982505555

For now I would suggest adding an ellipsis to the doctest to obscure the complex part (e.g write the output as 0.0145459669825000...). I've modified the patch to take this into account.

As everything else works fine, I would change the status to positive review but I want to hear your opinion about the modification first.

Changed 12 years ago by jvkersch

comment:5 Changed 12 years ago by jdemeyer

  • Authors set to Jeroen Demeyer
  • Reviewers set to Joris Vankerschaver

Joris, I agree with your change. I'm now testing on a 32-bit system...

comment:6 Changed 12 years ago by jdemeyer

  • Status changed from needs_review to positive_review

Tests successful, positive review as suggested by Joris.

comment:7 Changed 12 years ago by jvkersch

Sounds good, make ptestlong reports no failures here as well.

comment:8 Changed 12 years ago by jdemeyer

  • Merged in set to sage-4.6.1.alpha2
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.