Opened 5 years ago

Closed 3 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) Commit: 362c02d3d8a94a259c50d543281faadbb48631f7
Dependencies: Stopgaps:

Description (last modified by rws)

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 5 years ago by rws

  • Branch set to u/rws/special_values_of_bessel_functions

comment:2 Changed 5 years ago by rws

  • Authors set to Ralf Stephan
  • Cc chapoton added
  • Commit set to 191ca289d3c26813e0244251a983a62f2057347a
  • Description modified (diff)
  • Status changed from new to needs_review

New commits:

191ca2817678: special values of Bessel functions

comment:3 follow-up: Changed 5 years ago by kcrisman

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 5 years ago by rws

Replying to kcrisman:

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?

It does.

In #17130 Jeroen automatized handling of type in BuiltinFunction.

comment:5 Changed 5 years ago by kcrisman

I guess I didn't internalize that it would fix all future problems with this, very nice.

comment:6 follow-up: Changed 5 years ago by kcrisman

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: Changed 5 years ago by rws

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 for bessel_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: Changed 5 years ago by kcrisman

  • 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 5 years ago by rws

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 5 years ago by git

  • Commit changed from 191ca289d3c26813e0244251a983a62f2057347a to 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33

Branch pushed to git repo; I updated commit sha1. New commits:

c729ed9Merge branch 'develop' into t/17678/special_values_of_bessel_functions
1d27cb117678: return unsigned_infinity as special value; additions and fixes

comment:11 Changed 5 years ago by rws

  • 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 5 years ago by git

  • Commit changed from 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33 to c5f9994a07f76ece8a7b3ae2ed2a7df9fdcbbce6

Branch pushed to git repo; I updated commit sha1. New commits:

17d1aa117777: coerce unsigned infinity into SR
8025cd0Merge branch 'u/rws/unsigned_infinity_cannot_be_coerced_into_sr' of trac.sagemath.org:sage into t/17678/special_values_of_bessel_functions
c5f999417678: last fix

comment:13 Changed 5 years ago by rws

  • Status changed from needs_review to needs_work

There is a doctest fail in french_book that is not from #17777.

comment:14 Changed 5 years ago by git

  • Commit changed from c5f9994a07f76ece8a7b3ae2ed2a7df9fdcbbce6 to 882731883aa71630c6915cae75eafbc4777ec078

Branch pushed to git repo; I updated commit sha1. New commits:

dd4373eMerge branch 'develop' into t/17678/special_values_of_bessel_functions
f395febMerge branch 'develop' into t/17678/special_values_of_bessel_functions
0d05dd2Merge branch 'develop' into t/17777/unsigned_infinity_cannot_be_coerced_into_sr
3330a2fMerge branch 'develop' into t/17777/unsigned_infinity_cannot_be_coerced_into_sr
5609e5117777: fix doctests
9e0216eMerge branch 'u/rws/unsigned_infinity_cannot_be_coerced_into_sr' of trac.sagemath.org:sage into t/17678/special_values_of_bessel_functions
882731817678: doctest fixes

comment:15 Changed 5 years ago by rws

  • Branch changed from u/rws/special_values_of_bessel_functions to u/rws/17678

comment:16 Changed 5 years ago by rws

  • 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

Squashed it all into one commit.


New commits:

a215d3b17678: special values of Bessel functions

comment:17 in reply to: ↑ 7 ; follow-up: Changed 5 years ago by 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).

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: Changed 5 years ago by kcrisman

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 5 years ago by kcrisman

  • Reviewers set to Karl-Dieter Crisman
  • Status changed from needs_review to needs_info

comment:20 Changed 5 years ago by rws

  • Status changed from needs_info to needs_work

comment:21 in reply to: ↑ 18 Changed 5 years ago by rws

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 4 years ago by rws

  • Branch changed from u/rws/17678 to u/rws/17678-1

comment:23 Changed 4 years ago by rws

  • Commit changed from a215d3bc25cf6a9008e9807ed71f79d663c00236 to 8267479effcda361511e0b3c50dc2e86bbe94a54
  • Milestone changed from sage-6.6 to sage-6.10
  • Status changed from needs_work to needs_review

The infinity inconsistency is handled by #19439. Please review.


New commits:

8267479Merge branch 'u/rws/17678' of trac.sagemath.org:sage into tmp11

comment:24 in reply to: ↑ 18 ; follow-up: Changed 4 years ago by kcrisman

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 4 years ago by rws

  • Status changed from needs_review to needs_work

comment:26 Changed 3 years ago by arminstraub

  • Branch changed from u/rws/17678-1 to u/arminstraub/17678-1

comment:27 Changed 3 years ago by arminstraub

  • Authors changed from Ralf Stephan to Ralf Stephan, Armin Straub
  • 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 replaced n > 0 with n.real() > 0. I also changed the behaviour to only return unsigned_infinity if the real part of n is negative. The reason for these changes is that bessel_J(i, 0) should not be evaluated as 0, as it previously was, nor should it be evaluated as infinite. (Mathematica evaluates this value as Indeterminate.)
  • Likewise for bessel_I, for which the situation is equivalent.
  • Similarly, I changed Function_Bessel_Y._eval_ so that indeterminate values like bessel_Y(I,0) are not evaluated.
  • Likewise for bessel_K.

Doctesting the 5 involved files didn't reveal any troubles.


New commits:

02cb5f5Merge branch 'u/rws/17678-1' of git://trac.sagemath.org/sage into test_avoid_recompilation
bdd972f17678: adjust values at 0

comment:28 in reply to: ↑ 24 Changed 3 years ago by arminstraub

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: Changed 3 years ago by arminstraub

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: Changed 3 years ago by rws

  • 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 3 years ago by git

  • Commit changed from bdd972f6df2379729bc6caef0102974ef197a52c to 376ac2e646c7d27cc2bd10f3833a8a3eeb94a818

Branch pushed to git repo; I updated commit sha1. New commits:

376ac2e17678: fixed one doctest; avoid calling maxima for testing x==0

comment:32 in reply to: ↑ 30 ; follow-up: Changed 3 years ago by arminstraub

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 3 years ago by arminstraub

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: Changed 3 years ago by rws

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 for x because the index n 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 3 years ago by rws

  • Branch changed from u/arminstraub/17678-1 to public/17678

comment:36 follow-up: Changed 3 years ago by rws

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

362c02d17678: fix spurious doctest

comment:37 in reply to: ↑ 36 Changed 3 years ago by arminstraub

  • 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 with sage -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 3 years ago by arminstraub

Replying to rws:

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

Posted as #21201.

comment:39 Changed 3 years ago by vbraun

  • Branch changed from public/17678 to 362c02d3d8a94a259c50d543281faadbb48631f7
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.