Opened 11 years ago

Closed 11 years ago

# Bessel functions of real argument have small imaginary component when scipy is used

Reported by: Owned by: jvkersch jason, jkantor major sage-4.6.2 numerical bessel, scipy, complex jdemeyer, jason, jkantor sage-4.6.2.alpha3 Joris Vankerschaver Simon Spicer N/A

### Description

When the scipy algorithm is used to compute the Bessel I, J, Y, K function, the return value often has a small imaginary component (even though the argument is real):

```  sage: bessel_J(5, 1.5)
0.00179942176736061
sage: bessel_J(5, 1.5, algorithm='scipy')
0.00179942176736000 + 4.40731221543000e-19*I
```

This does not seem to be a problem with scipy, but is a result of the way in which scipy is called: the argument is first converted into a complex number, and then the corresponding function is called:

```  sage: import scipy.special
sage: scipy.special.jv(5, 1.5)
0.0017994217673606115
sage: scipy.special.jv(5, complex(1.5, 0.0))
(0.0017994217673606126+4.4073122154284958e-19j)
```

One annoying consequence of this behavior is that straightforward plotting of the bessel functions becomes problematic when the scipy algorithm is used (the command below works fine when another algorithm is used):

```  sage: Bessel(5, typ='J', algorithm='scipy').plot(1, 10)
verbose 0 (3989: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points.
verbose 0 (3989: plot.py, generate_plot_points) Last error message: 'unable to simplify to float approximation'
```

The present patch simply remedies this by adding a few lines to `bessel_{I, J, K, Y}` to check whether the input was real, and if so, to take the real part of the return value. This is definitely the easiest solution, but maybe a better one exists based on how Sage calls scipy.

### comment:1 Changed 11 years ago by jvkersch

• Status changed from new to needs_review

Updated doctests

### comment:2 Changed 11 years ago by spice

• Authors set to Joris Vankerschaver
• Reviewers set to Simon Spicer

This patch modifies three functions, but a new doctest was only added for one. This updated patch includes doctests for the remaining two. Otherwise, all checks out.

### comment:3 Changed 11 years ago by jvkersch

Hi Simon,

Thanks for spotting this omission and for uploading a new patch! I ran `make ptestlong` on it and all tests passed correctly, and the modifications look good to me too. So if you're willing to give the patch a positive review, I think you can go ahead and do so. Thanks again for your time and effort.

J.

### comment:4 follow-up: ↓ 5 Changed 11 years ago by spice

• Status changed from needs_review to positive_review

Hi Joris

Roger. I'm not sure if I'm allowed to positive review my own patch submission, but I've gone ahead anyway. 'Twas a pleasure reviewing your work :-)

Simon

### comment:5 in reply to: ↑ 4 Changed 11 years ago by jdemeyer

Hi Joris

Roger. I'm not sure if I'm allowed to positive review my own patch submission, but I've gone ahead anyway.

Well, Joris said it was okay.

### comment:6 follow-up: ↓ 7 Changed 11 years ago by jvkersch

I saw it more as you reviewing my contributions and me reviewing your changes, but now that we have the blessing of the release manager himself we're definitely in the clear :) Thanks guys.

### comment:7 in reply to: ↑ 6 Changed 11 years ago by jdemeyer

• Milestone set to sage-4.6.2

I saw it more as you reviewing my contributions and me reviewing your changes, but now that we have the blessing of the release manager himself we're definitely in the clear :)

Did I mention that these blessings come at a price? :-)

### comment:8 Changed 11 years ago by jdemeyer

• Status changed from positive_review to needs_work
• Work issues set to rebase

Needs to be rebased to sage-4.6.2.alpha2:

```patching file sage/functions/special.py
Hunk #5 FAILED at 875
Hunk #7 succeeded at 960 with fuzz 2 (offset 12 lines).
1 out of 7 hunks FAILED -- saving rejects to file sage/functions/special.py.rej
```

### comment:9 Changed 11 years ago by jvkersch

• Status changed from needs_work to needs_review

I've updated the patch. The file `sage/functions/special.py` passes doctests, and I'm currently running `make ptestlong`. I'm updating the status to needs_review since I don't want to give a positive review myself, but this should just be a formality.

Jeroen, why don't we have an Orval at one of the upcoming Sage days ;)

### comment:10 Changed 11 years ago by jdemeyer

• Work issues rebase deleted

### comment:11 Changed 11 years ago by jdemeyer

• Merged in set to sage-4.6.2.alpha3
• Resolution set to fixed
• Status changed from needs_review to closed

Looks good to me. Positive review.

Note: See TracTickets for help on using tickets.