make bessel_J symbolic

[jicama]

The motivation for this is

sage: plot(bessel_J(1, x), (x, 1, 10))
Traceback (most recent call last):
...
TypeError: Unable to convert x (='1-1/8*x^2+1/192*x^4-1/9216*x^6+1/737280*x^8-1/88473600*x^10+1/14863564800*x^12-1/3329438515200*x^14+1/958878292377600*x^16+O(x^17)') to real number.

The problem is that special functions, or at least bessel_J, can't currently be partially evaluated--that is, called with a SymbolicExpression as an argument. The model of good behavior is polylog, for which the above method produces a perfectly nice plot

sage: plot(polylog(1,x),(x,.1,.9)) #makes a fine plot

See discussion at http://groups.google.com/group/sage-support/browse_thread/thread/1b985b080ba2369e/7dee9eed953857f5#7dee9eed953857f5

---

Release manager:

Apply:

• attachment: trac_symbolic_bessel_v7.2.patch,
• attachment: trac_symbolic_bessel_v7-doctests.patch,
• attachment: bessel_2.2.patch
• attachment: trac_4102-bessel_doctest_fixes2.patch

Depends on #9880

CC: @jasongrout @kcrisman @benjaminfjones @sagetrac-dsm @burcin @eviatarbach

Component: calculus

Keywords: sd48

Author: Benjamin Jones, Eviatar Bach, Volker Braun

Reviewer: Karl-Dieter Crisman, Burcin Erocal

Merged: sage-5.11.rc0

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

[jasongrout]

Changed keywords from none to jason

[jasongrout]

Changed keywords from jason to none

[burcin]

comment:5

This ticket description was too broad. We have lots of tickets on fixing symbolic functions, see symbolics/functions for an overview.

See #3426 and #4230 for other tickets related to bessel functions. It would make sense to fix all these together, with a base class that handles the common properties of bessel functions (differentiation), and subclasses that implement numerical evaluation, etc. for each type.

[kcrisman]

comment:7

See also #10636, which I somehow never saw before. This sage-support thread yields an interesting related error:

sage: var('k')
k
sage: Z = sum( ((-1)^k*(x^(2*k+1))/factorial(2*k+1)),k,0,oo)
sage: Z
1/2*sqrt(pi)*sqrt(2)*sqrt(x)*bessel_j(1/2, x)
sage: Z(x=3)
1/2*sqrt(pi)*sqrt(2)*sqrt(3)*bessel_j(1/2, 3)
sage: Z(x=3).n()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:20822)()

TypeError: cannot evaluate symbolic expression numerically
sage: besse
bessel_I  bessel_J  bessel_K  bessel_Y  

Note that we apparently aren't yet converting Maxima's Bessel properly to 'our' uppercase version because of this.

[benjaminfjones]

comment:8

This wouldn't be hard to implement using #11143 as a template, but what to do about names? I guess new names for the symbolic Bessel functions should be chosen and deprecation notices added to the existing Bessel_J etc.

What new names do people like in the interim?

• bessel_J_symbolic
• bessel_J_sym
• bessel_Js

???

[burcin]

comment:9

Replying to @benjaminfjones:

What new names do people like in the interim?

• bessel_J_symbolic
• bessel_J_sym
• bessel_Js

???

What is wrong with taking over bessel_{I,J,K,Y}? Why do we need interim names?

Thanks for looking into this.

[kcrisman]

comment:10

What is wrong with taking over bessel_{I,J,K,Y}? Why do we need interim names?

I concur.

• bessel_J_symbolic

One could use that as the "symbolic" class one and then do the usual foo = Foo_Class()?

Thanks for looking into this.

Agreed!

[benjaminfjones]

Author: Benjamin Jones

[benjaminfjones]

comment:12

I thought I'd post some work in progress for all the Bessel function fans in the crowd. There is still work to be done, e.g.

• getting symbolic integration via Maxima to work
• adding maxima, scipy, PARI, etc. to the implemented backend systems
• more documentation and doctests

but at this point the patch applies cleanly to 5.4.1 and all the doctests in sage/functions/special.py pass.

If you are interested, have a look at the doctests included with the Bessel function for examples and let me know if you see any serious problems (other than the deficits I've listed above).

Thanks!

[benjaminfjones]

work in progress making Bessel functions symbolic

[kcrisman]

comment:13

Attachment: trac_symbolic_bessel.metaclass.2.patch.gz

Nice! A few comments of the type you solicited:

• Why typ and not type? Some Python reserved thing, maybe? But it looks like a typ-o to the (quickly reading) end user.
• I'd like to be able to plot f(x) = Bessel(0) but maybe that doesn't make sense? I guess a variable is necessary... anyway, just throwing it out there.
• f = maxima(Bessel(typ='K')(x,y)) turns out great, but does it convert back? Like f.derivative('y') is -(bessel_k(x+1,y)+bessel_k(x-1,y))/2, but does it then (when put in the _sage_ method) go back to "our" uppercase Bessel functions?
• Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.
• At least some of the error messages should be in doctests, maybe the ones with the wrong type and a non-implemented system.
• class_attrs['_conversions'] = {} --- what do we do with this in Maxima, then? Maybe it's better to raise an error; Maxima tends to otherwise just take things as new variables, which could be dangerous.
• How many of the currently-deleted doctests do you think would be worth preserving in the long run? Any deprecation needed here?
  Anyway, clearly a lot of planning and looks very promising!

[burcin]

comment:14

Thanks for the patch. Bessel functions are symbolic with so few lines of code. Amazing. :)

I like the metaclass idea. That never occured to me before as an option for functions with parameters, like the order here. I have a few questions:

• what is the advantage of creating a new symbolic function for each (type, order) pair, instead of having a generic function for each type that takes the order as a parameter? The latter would be similar to the design of the psi() function in GiNaC/Pynac.

• wouldn't this approach make life harder if in the future we want to add exact evaluation method for Bessel_J only?

I am also concerned about blowing up the symbolic function registry with these instances. Note that we keep an entry in a C++ list with a pointer to all the possible custom methods, and a Python dictionary with a function -> Pynac id mapping for each subclass of BuiltinFunction that is instantiated.

BTW, I suggest moving this code to a new file sage/functions/bessel.py.

[benjaminfjones]

comment:15

Replying to @kcrisman:

Thanks for the comments. I finally found some time to get back to this ticket :)

Nice! A few comments of the type you solicited:

• Why typ and not type? Some Python reserved thing, maybe? But it looks like a typ-o to the (quickly reading) end user.

That was intentional, I was trying not to shadow the builtin name.

• I'd like to be able to plot f(x) = Bessel(0) but maybe that doesn't make sense? I guess a variable is necessary... anyway, just throwing it out there.

I think f(x) = Bessel(0,x) is better, I don't like x being implicit on the right hand side.

• f = maxima(Bessel(typ='K')(x,y)) turns out great, but does it convert back? Like f.derivative('y') is -(bessel_k(x+1,y)+bessel_k(x-1,y))/2, but does it then (when put in the _sage_ method) go back to "our" uppercase Bessel functions?

Good question, I haven't tested it out. I think I'm going to rewrite most of the code anyway so I'll keep this in mind.

• Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.

???

• At least some of the error messages should be in doctests, maybe the ones with the wrong type and a non-implemented system.

Good point, I'll make sure to doctest the exceptions.

• class_attrs['_conversions'] = {} --- what do we do with this in Maxima, then? Maybe it's better to raise an error; Maxima tends to otherwise just take things as new variables, which could be dangerous.

I don't know about this. I had this in the case of the single parameter functions like Bessel(1,'K') that Maxima doesn't have equivalents for (it just has the two parameter functions). I think I'll scrap the idea of having all infinitely many one-parameter Bessel functions registered as their own symbolic functions (see Burcin's comments below).

• How many of the currently-deleted doctests do you think would be worth preserving in the long run? Any deprecation needed here?

Ideally all of them. Some will need to be modified because of the change in numerical back-end.

Anyway, clearly a lot of planning and looks very promising!

[benjaminfjones]

comment:16

Replying to @burcin:

Thanks for looking at it!

Thanks for the patch. Bessel functions are symbolic with so few lines of code. Amazing. :)

I like the metaclass idea. That never occured to me before as an option for functions with parameters, like the order here. I have a few questions:

• what is the advantage of creating a new symbolic function for each (type, order) pair, instead of having a generic function for each type that takes the order as a parameter? The latter would be similar to the design of the psi() function in GiNaC/Pynac.

The advantage I had in mind was just code reuse. In hindsight though, this approach makes the code more complicated and less maintainable that it needs to be. I'm going to refactor it into four generic functions in sage/functions/bessel.py as you suggest.

• wouldn't this approach make life harder if in the future we want to add exact evaluation method for Bessel_J only?

Good point..

I am also concerned about blowing up the symbolic function registry with these instances. Note that we keep an entry in a C++ list with a pointer to all the possible custom methods, and a Python dictionary with a function -> Pynac id mapping for each subclass of BuiltinFunction that is instantiated.

Also a good point. This occurred to me, but I didn't think through the consequences very much. I can see it being a problem if a user can inadvertently create a very large number of symbolic functions at runtime by doing something innocuous like:

sage: point([ (k, Bessel(k, 'J')(pi)) for k in range(1000) ])
... (1000 new symbolic function classes created!)

BTW, I suggest moving this code to a new file sage/functions/bessel.py.

Good idea. New patch coming soon...

[kcrisman]

comment:17

• Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.

???

Such as this and this. Just for forward-compatibility (instead of the percent business). Problem is that when I looked for ways to get around braces "naturally" occurring in LaTeX, there weren't necessarily a lot of them. (Ways, that is.)

[benjaminfjones]

comment:18

Replying to @kcrisman:

• Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.

???

Such as this and this. Just for forward-compatibility (instead of the percent business). Problem is that when I looked for ways to get around braces "naturally" occurring in LaTeX, there weren't necessarily a lot of them. (Ways, that is.)

Oh, right. I just didn't see what in the code you were referring to. I see now.

Anyway, to one of your earlier questions, with the new code I'm about to post the following conversions to and from Maxima work great:

sage: mb = maxima(bessel_I(1,x))
sage: mb.sage()                 
bessel_I(1, x)
sage:  x,y = var('x,y')
sage:  f = maxima(Bessel(typ='K')(x,y))
sage: f.derivative('x')
%pi*csc(%pi*x)*('diff(bessel_i(-x,y),x,1)-'diff(bessel_i(x,y),x,1))/2-%pi*bessel_k(x,y)*cot(%pi*x)
sage: m = f.derivative('x')
sage: m.sage()
-1/2*(x*bessel_I(-x, y)/y - x*bessel_I(x, y)/y + bessel_I(-x - 1, y) + bessel_I(x - 1, y))*pi*csc(pi*x) - pi*cot(pi*x)*bessel_K(x, y)

[benjaminfjones]

Attachment: trac_symbolic_bessel_v2.patch.gz

more work in progress

[kcrisman]

comment:19

Just as an FYI, this ask.sagemath question wants lots and lots of precision on Bessel Y. Which I think will be provided here via mpmath - just pointing out we should doctest it.

[kcrisman]

comment:20

Also, I think that #3426 and #4230 probably would be solved by a successful resolution of this ticket. Let's make sure to include any suggested (and useful) doctests from there here.

[benjaminfjones]

comment:21

Yet more work in progress, added lots of doctests, in particular to show that problems in #4230 and #3426 are fixed by this patch. One feature this patch adds is the ability to solve
Bessel's diffeq (using Maxima) and get a usable symbolic solution back (see the doctests in the module level docstring and in the Bessel() function).

Added brief mathematical exposition on the module docstring as well as some usage EXAMPLES.

I think the patch is more or less ready for initial review. All tests withing the new file sage/functions/bessel.py pass on my machine and sage-5.6.beta3, but there are several doctests failing in other places that reference the Bessel functions. I'll work on patching those up for the next iteration.

[benjaminfjones]

more work in progress, v3

[kcrisman]

comment:22

Attachment: trac_symbolic_bessel_v3.patch.gz

Dumb comments since I don't have time for proper review - and more importantly, there are no obvious horrible things (though I haven't gone in depth with branches yet). If all this really works and there are no typos, I think this will be a VERY nice addition.

• unqiue typo
• Verfify typo
• increasin typo
• I don't know why, but the math following "It follows from Bessel's..." doesn't look right in the doc (the :math: directive is not parsed)
• Is \operatorname{bessel_I}(\alpha, x) standard or should we just just the I sub whatever?
• `test_relation()` needs to be in double back ticks or have a link to the appropriate place in the ref manual
• Trac tickets should use :trac: syntax
• Does f(x) = Bessel(0): plot(f, (x, 1, 10)) work, or must one use Bessel(0)(x)? The example after that makes it look like maybe not...
• Bessel Y and Bessel K need a little filling out - and one of the :math: directives doesn't show up at all, the other isn't parsed right for some reason again
• I'd personally get rid of the whitespace changes in sage/functions/special.py - unlikely to have an effect, but not really necessary.
• How should we deal with the removal of the "maxima" and "pari" arguments? I'm not sure if it's really feasible to have a deprecation period for that, but we should discuss it.
• I suppose we don't need to keep all old tests, but some of them are rather different and might be worth putting in a TESTS section somewhere, just so that we don't have some unexpected regression only they catch.
• The switch to alpha from nu - I would rather not deprecate this either, but in principle someone could have used it as a keyword argument in the past...

[kcrisman]

Reviewer: Karl-Dieter Crisman

[burcin]

comment:76

I uploaded a new patch to fix doctest failures caused by #9880.

[kcrisman]

comment:77

Looks good and passes tests.

[jdemeyer]

comment:78

On silius (ia64):

sage -t --long devel/sage/sage/functions/bessel.py
**********************************************************************
File "devel/sage/sage/functions/bessel.py", line 258, in sage.functions.bessel.Function_Bessel_J
Failed example:
    A[0]
Expected:
    0.44005058574493355
Got:
    0.44005058574493366
**********************************************************************

[eviatarbach]

comment:79

Attachment: trac_4102-bessel_doctest_fixes2.patch.gz

New patch, adding a tolerance for the integral which is higher than the maximum error given by GSL.

[burcin]

comment:80

Thanks!

[jdemeyer]

Merged: sage-5.11.rc0

[eviatarbach]

comment:82

This is wrong:

sage: var('nu z')
(nu, z)
sage: bessel_J(nu, z).diff(nu)
-1/2*bessel_J(nu + 1, z) + 1/2*bessel_J(nu - 1, z)
sage: bessel_J(nu, z).diff(z)
-1/2*bessel_J(nu + 1, z) + 1/2*bessel_J(nu - 1, z)

[eviatarbach]

Attachment: trac4102_diff.patch.gz

[eviatarbach]

comment:83

New patch fixes this issue. Can this be merged in 5.11?

[jdemeyer]

comment:84

Replying to @eviatarbach:

New patch fixes this issue. Can this be merged in 5.11?

No, you should make a follow-up ticket for this.

[eviatarbach]

comment:85

Oh okay, then can this ticket be removed from 5.11? I'm just worried about having mathematically incorrect answers in the release.

[jdemeyer]

comment:86

Replying to @eviatarbach:

New patch fixes this issue. Can this be merged in 5.11?

Sorry, I really meant: perhaps yes it can be in sage-5.11, but in any case there needs to be a follow-up ticket (make it milestone: sage-5.11 and priority: blocker) and the new patch needs to be reviewed.

[jdemeyer]

comment:87

Also, add a doctest for the NotImplementedError (like the tests from [comment:82])

[eviatarbach]

comment:88

Okay, thank you. This is now #15019.