Opened 12 years ago

Closed 12 years ago

# Bessel method forces bessel_Y to use pari algorithm

Reported by: Owned by: jvkersch jason, jkantor major sage-4.6.1 numerical functions, bessel, pari sage-4.6.1.alpha2 Jeroen Demeyer Joris Vankerschaver N/A

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

### 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`.

### 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.

### 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.