Opened 9 years ago
Closed 9 years ago
#10251 closed defect (fixed)
Bessel functions of real argument have small imaginary component when scipy is used
Reported by: | jvkersch | Owned by: | jason, jkantor |
---|---|---|---|
Priority: | major | Milestone: | sage-4.6.2 |
Component: | numerical | Keywords: | bessel, scipy, complex |
Cc: | jdemeyer, jason, jkantor | Merged in: | sage-4.6.2.alpha3 |
Authors: | Joris Vankerschaver | Reviewers: | Simon Spicer |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
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.
Attachments (3)
Change History (14)
Changed 9 years ago by
comment:1 Changed 9 years ago by
- Status changed from new to needs_review
Changed 9 years ago by
comment:2 Changed 9 years ago by
- 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 9 years ago by
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 9 years ago by
- 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 9 years ago by
Replying to spice:
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 9 years ago by
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 9 years ago by
- Milestone set to sage-4.6.2
Replying to 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 :)
Did I mention that these blessings come at a price? :-)
comment:8 Changed 9 years ago by
- 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
Changed 9 years ago by
comment:9 Changed 9 years ago by
- 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 9 years ago by
- Work issues rebase deleted
comment:11 Changed 9 years ago by
- 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.
Updated doctests