special values of Bessel functions

[rwst]

At the moment everything Bessel is sent to mpmath and returns floating point but I_n(0) and J_n(0) are in (0,1). Also all I/J/K/Y_n(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)

CC: @fchapoton @fredrik-johansson

Component: symbolics

Author: Ralf Stephan, Armin Straub

Branch/Commit: 362c02d

Reviewer: Karl-Dieter Crisman, Ralf Stephan, Armin Straub

Issue created by migration from https://trac.sagemath.org/ticket/17678

[rwst]

Branch: u/rws/special_values_of_bessel_functions

[rwst]

Author: Ralf Stephan

[rwst]

New commits:

191ca28

17678: special values of Bessel functions

[rwst]

Commit: 191ca28

[kcrisman]

comment:3

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

[rwst]

comment:4

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.

[kcrisman]

comment:5

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

[kcrisman]

comment:6

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.

[rwst]

comment:7

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 Y_n(z)=(J_n(z)cos(npi)-J_-n(z))/sin(n*pi) (Abramowitz and Stegun 1972, p. 358).

[kcrisman]

comment:8

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?

[rwst]

comment:9

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.

[sagetrac-git]

Changed commit from 191ca28 to 1d27cb1

[sagetrac-git]

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

c729ed9	Merge branch 'develop' into t/17678/special_values_of_bessel_functions
1d27cb1	17678: return unsigned_infinity as special value; additions and fixes

[rwst]

Dependencies: #17777

[rwst]

comment:11

This uncovered a bug so we depend on #17777 as soon as it is resolved.

[sagetrac-git]

Changed commit from 1d27cb1 to c5f9994

[sagetrac-git]

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

17d1aa1	17777: coerce unsigned infinity into SR
8025cd0	Merge branch 'u/rws/unsigned_infinity_cannot_be_coerced_into_sr' of trac.sagemath.org:sage into t/17678/special_values_of_bessel_functions
c5f9994	17678: last fix

[rwst]

comment:13

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

[arminstraub]

Changed branch from u/rws/17678-1 to u/arminstraub/17678-1

[arminstraub]

comment:27

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:

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

[arminstraub]

Changed author from Ralf Stephan to Ralf Stephan, Armin Straub

[arminstraub]

Changed commit from 8267479 to bdd972f

[arminstraub]

comment:28

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).

[arminstraub]

comment:29

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?

[rwst]

Changed reviewer from Karl-Dieter Crisman to Karl-Dieter Crisman, Ralf Stephan

[rwst]

Changed dependencies from #17777 to none

[rwst]

comment:30

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.

[sagetrac-git]

Changed commit from bdd972f to 376ac2e

[sagetrac-git]

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

376ac2e

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

[arminstraub]

comment:32

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!

[arminstraub]

comment:33

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?

[rwst]

comment:34

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

[rwst]

Changed branch from u/arminstraub/17678-1 to public/17678

[rwst]

comment:36

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

[rwst]

Changed commit from 376ac2e to 362c02d

[arminstraub]

Changed reviewer from Karl-Dieter Crisman, Ralf Stephan to Karl-Dieter Crisman, Ralf Stephan, Armin Straub

[arminstraub]

comment:37

Looks good! Set to positive review.

Replying to @rwst:

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.

[arminstraub]

comment:38

Replying to @rwst:

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.

[vbraun]

Changed branch from public/17678 to 362c02d