Opened 8 years ago

Closed 6 years ago

# special values of Bessel functions

Reported by: Owned by: rws major sage-7.4 symbolics chapoton, fredrik.johansson Ralf Stephan, Armin Straub Karl-Dieter Crisman, Ralf Stephan, Armin Straub N/A 362c02d 362c02d3d8a94a259c50d543281faadbb48631f7

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

### comment:1 Changed 8 years ago by rws

• Branch set to u/rws/special_values_of_bessel_functions

### comment:2 Changed 8 years ago by rws

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

New commits:

 ​191ca28 `17678: special values of Bessel functions`

### comment:3 follow-up: ↓ 4 Changed 8 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 8 years ago by rws

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 8 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: ↓ 7 Changed 8 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: ↓ 17 Changed 8 years ago by rws

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

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

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 git

• Commit changed from 191ca289d3c26813e0244251a983a62f2057347a to 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33

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`

### comment:11 Changed 7 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 7 years ago by git

• Commit changed from 1d27cb1e1d45289b2d27a9e7b3b49ef1cda3ba33 to c5f9994a07f76ece8a7b3ae2ed2a7df9fdcbbce6

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`

### comment:13 Changed 7 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 7 years ago by git

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

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

### comment:16 Changed 7 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:

 ​a215d3b `17678: special values of Bessel functions`

### comment:17 in reply to: ↑ 7 ; follow-up: ↓ 18 Changed 7 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: ↓ 21 ↓ 24 Changed 7 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).

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 kcrisman

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

### comment:20 Changed 7 years ago by rws

• Status changed from needs_info to needs_work

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

```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 rws

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

### comment:23 Changed 7 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:

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

### comment:24 in reply to: ↑ 18 ; follow-up: ↓ 28 Changed 7 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 6 years ago by rws

• Status changed from needs_review to needs_work

### comment:26 Changed 6 years ago by arminstraub

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

### comment:27 Changed 6 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:

 ​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 arminstraub

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 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: ↓ 32 Changed 6 years ago by rws

• Dependencies #17777 deleted
• Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Ralf Stephan

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 git

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

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 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: ↓ 38 Changed 6 years ago by rws

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

### comment:35 Changed 6 years ago by rws

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

### comment:36 follow-up: ↓ 37 Changed 6 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:

 ​362c02d `17678: fix spurious doctest`

### comment:37 in reply to: ↑ 36 Changed 6 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.

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

On the other hand, wouldn't it be convenient to have a global function `is_trivial_zero`...