special values of Bessel functions

[rws]

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)

[rws]

New commits:

​191ca28

17678: special values of Bessel functions

[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

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

[kcrisman]

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

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

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

[kcrisman]

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?

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

[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

[rws]

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

[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

[rws]

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

[git]

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

[rws]

Squashed it all into one commit.

---

New commits:

​a215d3b

17678: special values of Bessel functions

[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 Y_n(z)=(J_n(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.

[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 Y_n(z)=(J_n(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?

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

[rws]

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

---

New commits:

​8267479

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

[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)
+        """

[arminstraub]

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]

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]

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?

[rws]

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.

[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]

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]

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?

[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

[rws]

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

[arminstraub]

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.

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