Opened 8 years ago
Closed 6 years ago
#17678 closed enhancement (fixed)
special values of Bessel functions
Reported by: | rws | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-7.4 |
Component: | symbolics | Keywords: | |
Cc: | chapoton, fredrik.johansson | Merged in: | |
Authors: | Ralf Stephan, Armin Straub | Reviewers: | Karl-Dieter Crisman, Ralf Stephan, Armin Straub |
Report Upstream: | N/A | Work issues: | |
Branch: | 362c02d (Commits, GitHub, GitLab) | Commit: | 362c02d3d8a94a259c50d543281faadbb48631f7 |
Dependencies: | Stopgaps: |
Description (last modified by )
At the moment everything Bessel is sent to mpmath and returns floating point but In(0) and Jn(0) are in (0,1). Also all I/J/K/Yn(x) with n in (-1/2,1/2) is a simple expression.
For the latter see p.10 of http://www.math.psu.edu/papikian/Kreh.pdf
This allows exact computation of series:
sage: bessel_I(2,x).series(x,10) 1/8*x^2 + 1/96*x^4 + 1/3072*x^6 + 1/184320*x^8 + Order(x^10)
Change History (39)
comment:1 Changed 8 years ago by
- Branch set to u/rws/special_values_of_bessel_functions
comment:2 Changed 8 years ago by
- Cc chapoton added
- Commit set to 191ca289d3c26813e0244251a983a62f2057347a
- Description modified (diff)
- Status changed from new to needs_review
comment:3 follow-up: ↓ 4 Changed 8 years ago by
I like the general idea of this ticket! What happens in the following scenario
+ if x == 0: + if n == 0: + return ZZ(1)
for
sage: bessel_J(0.0, 0.0) 1
should it instead return 1.0
? I feel like we have a lot of code in Sage doing type conversions of this kind smartly - sadly, I was not the one gifted to write it. Maybe something like saving the parents of the input and if one is inexact then so is the output, something along those lines... we had some issues with sin(0)
versus sin(0.0)
even, I think; currently
sage: sin(0) 0 sage: sin(0.0) 0.000000000000000 sage: sin(float(0)) 0.0 sage: sin(RDF(0)) 0.0 sage: sin(complex(0)) 0j
comment:4 in reply to: ↑ 3 Changed 8 years ago by
comment:5 Changed 8 years ago by
I guess I didn't internalize that it would fix all future problems with this, very nice.
comment:6 follow-up: ↓ 7 Changed 8 years ago by
These identities are also available (or derivable) from Wikipedia and Mathworld, so we are in good shape.
Question: W|A claims that one has the complex infinity, not positive infinity, for some of the negative ones like bessel_J(-5/2, 0)
or for bessel_I
. I don't know what to make of that, though. Also, does bessel_Y
have an analogous special value for x=0, negative n?
Otherwise looks good.
comment:7 in reply to: ↑ 6 ; follow-up: ↓ 17 Changed 8 years ago by
Replying to kcrisman:
Question: W|A claims that one has the complex infinity, not positive infinity, for some of the negative ones like
bessel_J(-5/2, 0)
or forbessel_I
. I don't know what to make of that, though.
I got those values from mpmath and just tried to post about that to the mpmath group, but not yet approved.
Also, does
bessel_Y
have an analogous special value for x=0, negative n?
Ah, I missed that. It should be easily derived from Bessel_J
with Yn(z)=(Jn(z)*cos(n*pi)-J-n(z))/sin(n*pi) (Abramowitz and Stegun 1972, p. 358).
comment:8 follow-up: ↓ 9 Changed 7 years ago by
- Cc fredrik.johansson added
- Status changed from needs_review to needs_work
Ah, I missed that. It should be easily derived from
Yes, I figured - but should it be included? Sorry for not being clear.
Did you hear back from Fredrik/mpmath?
comment:9 in reply to: ↑ 8 Changed 7 years ago by
Replying to kcrisman:
Ah, I missed that. It should be easily derived from
Yes, I figured - but should it be included?
Of course!
Did you hear back from Fredrik/mpmath?
https://groups.google.com/forum/?hl=en#!topic/mpmath/FJqtBMNhYFo So, IMO the mpmath behaviour fits its numerical requirements but we should use "zoo" because it's more symbolically correct, and because we have it.
comment:10 Changed 7 years ago by
- Commit changed from 191ca289d3c26813e0244251a983a62f2057347a to 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33
comment:11 Changed 7 years ago by
- Dependencies set to #17777
- Status changed from needs_work to needs_review
This uncovered a bug so we depend on #17777 as soon as it is resolved.
comment:12 Changed 7 years ago by
- Commit changed from 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33 to c5f9994a07f76ece8a7b3ae2ed2a7df9fdcbbce6
comment:13 Changed 7 years ago by
- Status changed from needs_review to needs_work
There is a doctest fail in french_book that is not from #17777.
comment:14 Changed 7 years ago by
- Commit changed from c5f9994a07f76ece8a7b3ae2ed2a7df9fdcbbce6 to 882731883aa71630c6915cae75eafbc4777ec078
Branch pushed to git repo; I updated commit sha1. New commits:
dd4373e | Merge branch 'develop' into t/17678/special_values_of_bessel_functions
|
f395feb | Merge branch 'develop' into t/17678/special_values_of_bessel_functions
|
0d05dd2 | Merge branch 'develop' into t/17777/unsigned_infinity_cannot_be_coerced_into_sr
|
3330a2f | Merge branch 'develop' into t/17777/unsigned_infinity_cannot_be_coerced_into_sr
|
5609e51 | 17777: fix doctests
|
9e0216e | Merge branch 'u/rws/unsigned_infinity_cannot_be_coerced_into_sr' of trac.sagemath.org:sage into t/17678/special_values_of_bessel_functions
|
8827318 | 17678: doctest fixes
|
comment:15 Changed 7 years ago by
- Branch changed from u/rws/special_values_of_bessel_functions to u/rws/17678
comment:16 Changed 7 years ago by
- Commit changed from 882731883aa71630c6915cae75eafbc4777ec078 to a215d3bc25cf6a9008e9807ed71f79d663c00236
- Milestone changed from sage-6.5 to sage-6.6
- Priority changed from minor to major
- Status changed from needs_work to needs_review
comment:17 in reply to: ↑ 7 ; follow-up: ↓ 18 Changed 7 years ago by
Also, does
bessel_Y
have an analogous special value for x=0, negative n?Ah, I missed that. It should be easily derived from
Bessel_J
with Yn(z)=(Jn(z)*cos(n*pi)-J-n(z))/sin(n*pi) (Abramowitz and Stegun 1972, p. 358).
I didn't see anything about this in the many commits here, haven't looked at the 'squashed' one but assume it's not there either. Otherwise this just requires (from my point of view) some testing.
comment:18 in reply to: ↑ 17 ; follow-ups: ↓ 21 ↓ 24 Changed 7 years ago by
Replying to kcrisman:
Also, does
bessel_Y
have an analogous special value for x=0, negative n?Ah, I missed that. It should be easily derived from
Bessel_J
with Yn(z)=(Jn(z)*cos(n*pi)-J-n(z))/sin(n*pi) (Abramowitz and Stegun 1972, p. 358).
Also:
sage: bessel_J(-5./2,0.) +infinity sage: type(_) <type 'sage.rings.real_mpfr.RealNumber'>
I'm not sure what the "right" resolution is here, since MPFR may not have an unsigned infinity.
Patchbot is happy and I don't expect my own last test run to go bad, so these two things are all that remains to resolve, perhaps they are even non-issues?
comment:19 Changed 7 years ago by
- Reviewers set to Karl-Dieter Crisman
- Status changed from needs_review to needs_info
comment:20 Changed 7 years ago by
- Status changed from needs_info to needs_work
comment:21 in reply to: ↑ 18 Changed 7 years ago by
Replying to kcrisman:
sage: bessel_J(-5./2,0.) +infinity sage: type(_) <type 'sage.rings.real_mpfr.RealNumber'>I'm not sure what the "right" resolution is here, since MPFR may not have an unsigned infinity.
This appears to be widespread, e.g. I just found polylog(1.,1.)
vs. polylog(1,1)
which is due to zeta(1.)
vs. zeta(1)
.
comment:22 Changed 7 years ago by
- Branch changed from u/rws/17678 to u/rws/17678-1
comment:23 Changed 7 years ago by
- Commit changed from a215d3bc25cf6a9008e9807ed71f79d663c00236 to 8267479effcda361511e0b3c50dc2e86bbe94a54
- Milestone changed from sage-6.6 to sage-6.10
- Status changed from needs_work to needs_review
comment:24 in reply to: ↑ 18 ; follow-up: ↓ 28 Changed 7 years ago by
Thanks for opening that ticket, that is fine for that issue.
Still a problem here:
sage: bessel_Y(-3,0) Infinity sage: bessel_Y(3,0) Infinity sage: bessel_Y(3.1,0) -infinity sage: bessel_Y(3.2,0) -infinity sage: bessel_Y(-3.2,0) +infinity sage: bessel_Y(-3.3,0) +infinity sage: bessel_Y(-33/10,0) Infinity
Extra var('n')
?
+ def _eval_(self, n, x): + """ + EXAMPLES:: + + sage: n = var('n') + sage: bessel_Y(1, 0) + Infinity + sage: bessel_Y(1/2, x) + -sqrt(2)*sqrt(1/(pi*x))*cos(x) + sage: bessel_Y(-1/2, x) + sqrt(2)*sqrt(1/(pi*x))*sin(x) + """
comment:25 Changed 6 years ago by
- Status changed from needs_review to needs_work
comment:26 Changed 6 years ago by
- Branch changed from u/rws/17678-1 to u/arminstraub/17678-1
comment:27 Changed 6 years ago by
- Commit changed from 8267479effcda361511e0b3c50dc2e86bbe94a54 to bdd972f6df2379729bc6caef0102974ef197a52c
- Milestone changed from sage-6.10 to sage-7.4
- Status changed from needs_work to needs_review
We stumpled across the issue that bessel_J(0,x).series(x,3)
didn't have the expected result during a summer school on special functions, in which I was advertising Sage. It would be nice to report back to the students that the next version of Sage has this fixed :)
This is my first attempt at using git and I haven't used the trac server in many years, so please let me know if I messed something up or didn't follow best practices. The branch I pushed is supposed to be a merge with the most recent version of Sage, with the following additional changes:
- Removed the extra
var('n')
as noticed by Karl-Dieter. - In
Function_Bessel_J._eval_
, I replacedn > 0
withn.real() > 0
. I also changed the behaviour to only returnunsigned_infinity
if the real part ofn
is negative. The reason for these changes is thatbessel_J(i, 0)
should not be evaluated as0
, as it previously was, nor should it be evaluated as infinite. (Mathematica evaluates this value asIndeterminate
.) - Likewise for
bessel_I
, for which the situation is equivalent. - Similarly, I changed
Function_Bessel_Y._eval_
so that indeterminate values likebessel_Y(I,0)
are not evaluated. - Likewise for
bessel_K
.
Doctesting the 5 involved files didn't reveal any troubles.
New commits:
02cb5f5 | Merge branch 'u/rws/17678-1' of git://trac.sagemath.org/sage into test_avoid_recompilation
|
bdd972f | 17678: adjust values at 0
|
comment:28 in reply to: ↑ 24 Changed 6 years ago by
Replying to kcrisman:
Still a problem here:
sage: bessel_Y(-3,0) Infinity sage: bessel_Y(3,0) Infinity sage: bessel_Y(3.1,0) -infinity sage: bessel_Y(3.2,0) -infinity sage: bessel_Y(-3.2,0) +infinity sage: bessel_Y(-3.3,0) +infinity sage: bessel_Y(-33/10,0) Infinity
Indeed, the behaviour for values of Bessel functions at zero is not consistent between symbolic and numeric input. The symbolic evaluations as provided by this ticket seem appropriate to me (and, as far as I have tested, also agree with the values that Mathematica produces).
On the other hand, the numerical evaluation of Bessel functions is currently outsourced to mpmath
. I don't know about the conventions that mpmath
is using when reporting the corresponding values for numeric x
(for instance, for bessel_Y(1/2,x) == -sqrt(2)*sqrt(1/(pi*x))*cos(x)
mpmath produces mpmath.bessely(0.5,0) == mpf('-inf')
, which results in bessel_Y(0.5,0) == -infinity
, which is clearly not the limiting value when x approaches zero on the imaginary axis). To achieve perfect consistency, one could modify the _evalf_
function to not let mpmath
handle these cases. This seems, however, not like an ideal solution.
I would suggest that resolving this inconsistency is better suited for a different ticket (or, if desired, changing the behaviour within mpmath
itself).
comment:29 follow-up: ↓ 30 Changed 6 years ago by
By the way, what is the preferred approach of Sage to the following?
When the index of the Bessel functions is a half-integer, they can be written in terms of elementary functions. This is currently implemented for indices 1/2 and -1/2 only. Would it be desirable to likewise implement the case of general half-integer indices? Or, would it be better to leave, say, bessel_J(21/2, x)
unevaluated and wait for the user to, somehow (how?), explicitly ask for a simplification?
comment:30 in reply to: ↑ 29 ; follow-up: ↓ 32 Changed 6 years ago by
- Dependencies #17777 deleted
- Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Ralf Stephan
Replying to arminstraub:
It would be nice to report back to the students that the next version of Sage has this fixed :)
Indeed, but the discussion focused a bit on the dependency on #17777, which I agree however does not exist, i.e. both issues are mutually independent.
Your additions look fine, consider them reviewed. There is one failing doctest due to simplification in src/sage/interfaces/maxima_lib.py
. A minor caveat is also in my previous branch, all the comparisons n,x == 0
should rather read n,x.is_trivial_zero()
because any nonnumerical symbolic expression will trigger the expensive proof machinery in Expression.__nonzero__()
. So if you could please change these two minor issues you can set positive yourself. Thanks for your work.
By the way, what is the preferred approach of Sage to the following?
When the index of the Bessel functions is a half-integer, they can be written in terms of elementary functions. This is currently implemented for indices 1/2 and -1/2 only. Would it be desirable to likewise implement the case of general half-integer indices? Or, would it be better to leave, say,
bessel_J(21/2, x)
unevaluated and wait for the user to, somehow (how?), explicitly ask for a simplification?
This was not formalized up to now, the general behaviour was to give such conversions if the result is both more elementary and not very complicated. Full implementation can be done in a Expression.expand_bessel()
function which would then, at some time later, be part of a general rewrite()
tool (#10137). So to your question, automatic expansion yes, as long as the output doesn't exceed---say one line, and in a dedicated expand function always welcome.
comment:31 Changed 6 years ago by
- Commit changed from bdd972f6df2379729bc6caef0102974ef197a52c to 376ac2e646c7d27cc2bd10f3833a8a3eeb94a818
Branch pushed to git repo; I updated commit sha1. New commits:
376ac2e | 17678: fixed one doctest; avoid calling maxima for testing x==0
|
comment:32 in reply to: ↑ 30 ; follow-up: ↓ 34 Changed 6 years ago by
Thanks for your swift and helpful feedback!
I have fixed the failing doctest (the entire example needed to be replaced).
You also mentioned to use sth.is_trivial_zero()
instead of sth == 0
, and the reasoning makes perfect sense. A minor trouble is that when calling, say, sage: bessel_J(5/2, 0)
, then x
and n
are not expressions and so don't have the property is_trivial_zero
. One workaround would be to use (isinstance(sth, Expression) and sth.is_trivial_zero()) or sth == 0
but that doesn't feel natural to me. Another workaround would be SR(sth).is_trivial_zero()
but that also doesn't feel right (and I couldn't find this being used in other Sage code). After some more looking, I found that the code for hypergeometric functions uses things like if not isinstance(z, Expression) and z == 0:
a couple of times. For now, that's what I used also (only for x
because the index n
shouldn't usually be a complicated expression, and because other functions distinguish similarly).
Since there was some choices to be made, I haven't set the ticket to positive review myself. Could you please have a quick look!
comment:33 Changed 6 years ago by
On the other hand, wouldn't it be convenient to have a global function is_trivial_zero
for the purpose of comparing with zero, without working too hard to simplify?
def is_trivial_zero(x): try: return x.is_trivial_zero() except AttributeError: return x == 0
One place for such a function could be sage/misc/functional.py
which already contains things like is_even
(which I used as a blueprint here).
Do you think it would be a good idea to have something like that? If so, should a ticket be opened for that?
comment:34 in reply to: ↑ 32 ; follow-up: ↓ 38 Changed 6 years ago by
Replying to arminstraub:
... to use
(isinstance(sth, Expression) and sth.is_trivial_zero()) or sth == 0
No, that calls __nonzero__
after is_trivial_zero()
returned False
. You need (isinstance(sth, Expression) and sth.is_trivial_zero()) or (not isinstance(sth, Expression) and sth == 0)
.
After some more looking, I found that the code for hypergeometric functions uses things like
if not isinstance(z, Expression) and z == 0:
a couple of times. For now, that's what I used also (only forx
because the indexn
shouldn't usually be a complicated expression, and because other functions distinguish similarly).
Well that misses the SR(0)
case and looks buggy.
On the other hand, wouldn't it be convenient to have a global function
is_trivial_zero
...
Try opening a ticket, but see also #17158
comment:35 Changed 6 years ago by
- Branch changed from u/arminstraub/17678-1 to public/17678
comment:36 follow-up: ↓ 37 Changed 6 years ago by
- Commit changed from 376ac2e646c7d27cc2bd10f3833a8a3eeb94a818 to 362c02d3d8a94a259c50d543281faadbb48631f7
Well that misses the
SR(0)
case and looks buggy.
No, that actually works. What's odd is that I get now an error in src/sage/tests/french_book/recequadiff.py
that only happens with sage -tp
. As it is spurious I changed the test. As your work is reviewed someone still needs to look at my commit. Could you please do that and then set positive?
New commits:
362c02d | 17678: fix spurious doctest
|
comment:37 in reply to: ↑ 36 Changed 6 years ago by
- Reviewers changed from Karl-Dieter Crisman, Ralf Stephan to Karl-Dieter Crisman, Ralf Stephan, Armin Straub
- Status changed from needs_review to positive_review
Looks good! Set to positive review.
Replying to rws:
Well that misses the
SR(0)
case and looks buggy.No, that actually works.
Yes, I was initially worried about that, too, when following the hypergeometric implementation.
What's odd is that I get now an error in
src/sage/tests/french_book/recequadiff.py
that only happens withsage -tp
. As it is spurious I changed the test.
Your solution seems fine. I am seeing the same spurious behavior when running the test versus running the code in a notebook.
You need
(isinstance(sth, Expression) and sth.is_trivial_zero()) or (not isinstance(sth, Expression) and sth == 0)
.
Oops, you are right! Definitely not a brief and convenient substitute for sth == 0
. I'll open a ticket suggesting a global is_trivial_zero
function.
comment:38 in reply to: ↑ 34 Changed 6 years ago by
comment:39 Changed 6 years ago by
- Branch changed from public/17678 to 362c02d3d8a94a259c50d543281faadbb48631f7
- Resolution set to fixed
- Status changed from positive_review to closed
New commits:
17678: special values of Bessel functions