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)

10251_bessel_scipy.patch (2.3 KB) - added by jvkersch 9 years ago.
10251_bessel_scipy.2.patch (3.3 KB) - added by spice 9 years ago.
Updated doctests
10251_rebase_bessel_scipy.2.patch (3.2 KB) - added by jvkersch 9 years ago.

Download all attachments as: .zip

Change History (14)

Changed 9 years ago by jvkersch

comment:1 Changed 9 years ago by jvkersch

  • Status changed from new to needs_review

Changed 9 years ago by spice

Updated doctests

comment:2 Changed 9 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 9 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: Changed 9 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 9 years ago by jdemeyer

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: Changed 9 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 9 years ago by jdemeyer

  • 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 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

Changed 9 years ago by jvkersch

comment:9 Changed 9 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 9 years ago by jdemeyer

  • Work issues rebase deleted

comment:11 Changed 9 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.