problems with Bessel functions

[Pablo De Napoli]

Hi

I'm trying to do some computations with Bessel functions using Sage.
Unfortunately, they don't seem to behave like other functions. For example:
to get the plot of the sine function over the interval (0,100)

plot(sin(x),(x,0,100))

works. However,

plot(bessel_J(0,x),(x,0,100))

does not. Gives the following error messages:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)

/home/pablo/sage/sage-4.0.rc0/<ipython console> in <module>()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/functions/special.pyc
in bessel_J(nu, z, algorithm, prec)
753 C = ComplexField(prec)
754 nu = C(nu)
--> 755 z = C(z)
756 K = C
757 if nu == 0:

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/rings/complex_field.pyc
in __call__(self, x, im)
256 if im is not None:
257 x = x, im
--> 258 return Parent.__call__(self, x)
259
260 def _element_constructor_(self, x):

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/parent.so
in sage.structure.parent.Parent.__call__
(sage/structure/parent.c:4121)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/coerce_maps.so
in sage.structure.coerce_maps.DefaultConvertMap_unique._call_
(sage/structure/coerce_maps.c:3064)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/coerce_maps.so
in sage.structure.coerce_maps._call_
(sage/structure/coerce_maps.c:2955)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/rings/complex_field.pyc
in _element_constructor_(self, x)
280
281 try:
--> 282 return x._complex_mpfr_field_( self )
283 except AttributeError:
284 pass

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so
in sage.symbolic.expression.Expression._complex_mpfr_field_
(sage/symbolic/expression.cpp:5573)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so
in sage.symbolic.expression.Expression._eval_self
(sage/symbolic/expression.cpp:5027)()

TypeError: Cannot evaluate symbolic expression to a numeric value.

In a similar way, you can for example compute the derivative of the
sine function by

sage: f(x)=sin(x)
sage: f.diff()
x |--> cos(x)

but if you try to define a function by

f(x)=bessel_J(0,x)

you get a similar error:

---------------------------------------------------------------------------
TypeError Traceback (most recent call last)

/home/pablo/sage/sage-4.0.rc0/<ipython console> in <module>()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/functions/special.pyc
in bessel_J(nu, z, algorithm, prec)
753 C = ComplexField(prec)
754 nu = C(nu)
--> 755 z = C(z)
756 K = C
757 if nu == 0:

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/rings/complex_field.pyc
in __call__(self, x, im)
256 if im is not None:
257 x = x, im
--> 258 return Parent.__call__(self, x)
259
260 def _element_constructor_(self, x):

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/parent.so
in sage.structure.parent.Parent.__call__
(sage/structure/parent.c:4121)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/coerce_maps.so
in sage.structure.coerce_maps.DefaultConvertMap_unique._call_
(sage/structure/coerce_maps.c:3064)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/structure/coerce_maps.so
in sage.structure.coerce_maps._call_
(sage/structure/coerce_maps.c:2955)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/rings/complex_field.pyc
in _element_constructor_(self, x)
280
281 try:
--> 282 return x._complex_mpfr_field_( self )
283 except AttributeError:
284 pass

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so
in sage.symbolic.expression.Expression._complex_mpfr_field_
(sage/symbolic/expression.cpp:5573)()

/home/pablo/sage/sage-4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so
in sage.symbolic.expression.Expression._eval_self
(sage/symbolic/expression.cpp:5027)()

TypeError: Cannot evaluate symbolic expression to a numeric value.

=========================================================

I think it would be important to provide some uniform interfase for all
functions.

regards
Pablo

[Dan Drake]

On Tue, 08 Dec 2009 at 10:44PM -0300, Pablo De Napoli wrote:
> I'm trying to do some computations with Bessel functions using Sage.
> Unfortunately, they don't seem to behave like other functions. For example:
> to get the plot of the sine function over the interval (0,100)
>
> plot(sin(x),(x,0,100))
>
> works. However,
>
> plot(bessel_J(0,x),(x,0,100))
>
> does not.

Part of the reason is that (I think) the Bessel functions are not
implemented as fully symbolic functions; they are more numerical in
nature under our current implementation.

An effective workaround for your problems is to make anonymous
functions:

plot(lambda x: bessel_J(0, x), (0, 100))

which gives plot() a regular Python function of one variable that it can
use.

Dan

--
--- Dan Drake
----- http://mathsci.kaist.ac.kr/~drake
-------

[Nick Alexander]

On 8-Dec-09, at 8:59 PM, Dan Drake wrote:

> On Tue, 08 Dec 2009 at 10:44PM -0300, Pablo De Napoli wrote:
>> I'm trying to do some computations with Bessel functions using Sage.
>> Unfortunately, they don't seem to behave like other functions. For
>> example:
>> to get the plot of the sine function over the interval (0,100)
>>
>> plot(sin(x),(x,0,100))
>>
>> works. However,
>>
>> plot(bessel_J(0,x),(x,0,100))
>>
>> does not.
>
> Part of the reason is that (I think) the Bessel functions are not
> implemented as fully symbolic functions; they are more numerical in
> nature under our current implementation.

Implementing a particular symbolic function is not outlandishly
difficult, thanks to the tireless work of Burcin Erocal and Mike
Hansen. (Apologies to any contributers I have forgotten.)

You need to subclass sage.symbolic.function.SFunction. I don't see
many examples, so here is a minimal one:

from sage.symbolic.function import SFunction
from sage.rings.all import RealField

class bessel_J_class(SFunction):
def __init__(self, *args, **kwds):
kwds['nargs'] = 2
kwds['evalf_func'] = self._evalf_func_
SFunction.__init__(self, "bessel_J", *args, **kwds)

def _evalf_func_(self, *args, **kwds):
prec = kwds['prec']
vals = [ arg.n(prec) for arg in args ]
v = bessel_J(*vals)
return RealField(prec)(v)

symbolic_bessel_J = bessel_J_class()

Then the following works for me:

sage: var('x')
sage: plot(symbolic_bessel_J(0, x), (x, 0, 100))

Nick

[Mike Hansen]

Hello,

On Wed, Dec 9, 2009 at 12:44 PM, Nick Alexander <ncale...@gmail.com> wrote:
> Implementing a particular symbolic function is not outlandishly
> difficult, thanks to the tireless work of Burcin Erocal and Mike
> Hansen.  (Apologies to any contributers I have forgotten.)
>
> You need to subclass sage.symbolic.function.SFunction.  I don't see
> many examples, so here is a minimal one:
>
> from sage.symbolic.function import SFunction
> from sage.rings.all import RealField
>
> class bessel_J_class(SFunction):
>     def __init__(self, *args, **kwds):
>         kwds['nargs'] = 2
>         kwds['evalf_func'] = self._evalf_func_
>         SFunction.__init__(self, "bessel_J", *args, **kwds)
>
>     def _evalf_func_(self, *args, **kwds):
>         prec = kwds['prec']
>         vals = [ arg.n(prec) for arg in args ]
>         v = bessel_J(*vals)
>         return RealField(prec)(v)
>
> symbolic_bessel_J = bessel_J_class()
>
> Then the following works for me:
>
> sage: var('x')
> sage: plot(symbolic_bessel_J(0, x), (x, 0, 100))

Things have changed a tiny bit in Sage 4.3 which should make adding
these types of functions a bit easier. See patch # 7490 for lots of
examples.

--Mike

[Robert Bradshaw]

Any thoughts on making a symbolic decorator, so one could do something
like

@sage.symbolic.function.symbolic
def my_func(x, n):
if x < 0: return 0
else: return exp(-1/x^n)

?

- Robert

[Burcin Erocal]

On Tue, 8 Dec 2009 21:57:30 -0800
Robert Bradshaw <robe...@math.washington.edu> wrote:

<snip>

> > On Wed, Dec 9, 2009 at 12:44 PM, Nick Alexander
> > <ncale...@gmail.com> wrote:

<snip>

> >> You need to subclass sage.symbolic.function.SFunction. I don't see
> >> many examples, so here is a minimal one:
> >>
> >> from sage.symbolic.function import SFunction
> >> from sage.rings.all import RealField
> >>
> >> class bessel_J_class(SFunction):
> >> def __init__(self, *args, **kwds):
> >> kwds['nargs'] = 2
> >> kwds['evalf_func'] = self._evalf_func_
> >> SFunction.__init__(self, "bessel_J", *args, **kwds)
> >>
> >> def _evalf_func_(self, *args, **kwds):
> >> prec = kwds['prec']
> >> vals = [ arg.n(prec) for arg in args ]
> >> v = bessel_J(*vals)
> >> return RealField(prec)(v)
> >>
> >> symbolic_bessel_J = bessel_J_class()
> >>
> >> Then the following works for me:
> >>
> >> sage: var('x')
> >> sage: plot(symbolic_bessel_J(0, x), (x, 0, 100))

This is also #1158 on trac:

http://trac.sagemath.org/sage_trac/ticket/1158

<snip>

> Any thoughts on making a symbolic decorator, so one could do
> something like
>
> @sage.symbolic.function.symbolic
> def my_func(x, n):
> if x < 0: return 0
> else: return exp(-1/x^n)
>
> ?

Great idea! I opened a ticket:

http://trac.sagemath.org/sage_trac/ticket/7636

I'll implement this first chance I get.


Cheers,
Burcin

[Jason Grout]

Robert Bradshaw wrote:

> @sage.symbolic.function.symbolic
> def my_func(x, n):
> if x < 0: return 0
> else: return exp(-1/x^n)

How about exposing the many other properties too:

@sage.symbolic.function.symbolic(latex_name="\phi",...)
def my_func(x,n):
...


which translates to:

def my_func_eval(x,n):

if x < 0: return 0

else ...

my_func = function(..., latex_name="\phi",...(any other options),
eval_f=my_func_eval)


Thanks,

Jason


--
Jason Grout

[Jason Grout]

Nick Alexander wrote:
> Implementing a particular symbolic function is not outlandishly
> difficult, thanks to the tireless work of Burcin Erocal and Mike
> Hansen. (Apologies to any contributers I have forgotten.)
>
> You need to subclass sage.symbolic.function.SFunction. I don't see
> many examples, so here is a minimal one:
>
> from sage.symbolic.function import SFunction
> from sage.rings.all import RealField
>
> class bessel_J_class(SFunction):
> def __init__(self, *args, **kwds):
> kwds['nargs'] = 2
> kwds['evalf_func'] = self._evalf_func_
> SFunction.__init__(self, "bessel_J", *args, **kwds)
>
> def _evalf_func_(self, *args, **kwds):
> prec = kwds['prec']
> vals = [ arg.n(prec) for arg in args ]
> v = bessel_J(*vals)
> return RealField(prec)(v)
>
> symbolic_bessel_J = bessel_J_class()
>
> Then the following works for me:
>
> sage: var('x')
> sage: plot(symbolic_bessel_J(0, x), (x, 0, 100))

Is there any difference between declaring the class above and doing the
following?

symbolic_bessel_J=function('symbolic_bessel_J',nargs=2,evalf_func=bessel_J)

This thread also might be useful to understand evalf:

http://groups.google.com/group/sage-devel/browse_frm/thread/983cd9ca5977dd2d/0c47026ce3914157?lnk=gst&q=partially#0c47026ce3914157

[Nathann Cohen]

Would it be possible, using this, to define a symbolic Sum ? Something
like Sum(Set([1,2,3,4,5,6]))...This would be extremely useful in
LP !!!

Nathann

> http://groups.google.com/group/sage-devel/browse_frm/thread/983cd9ca5...

[Mike Hansen]

On Wed, Dec 9, 2009 at 7:01 PM, Nathann Cohen <nathan...@gmail.com> wrote:
> Would it be possible, using this, to define a symbolic Sum ? Something
> like Sum(Set([1,2,3,4,5,6]))...This would be extremely useful in

What would the output of that be?

There is http://trac.sagemath.org/sage_trac/ticket/3587 , but I don't
know that's what you're after.

--Mike

[Pablo De Napoli]

Many thanks to everybody for your help.

some questions/remarks:

1) ¿Does every function needs to have two versions: a symbolic one and
a numerical
one?

2) the current implementation treats f(x)=sin(x) as a symbolic expression

sage: f(x)=sin(x)
sage: f
x |--> sin(x)
sage: type(f)
<type 'sage.symbolic.expression.Expression'>

However, mathematically a function is something different: for instance it
has a domain and a range (say: this function takes one real number as
argument, and returns another real number).

¿is there some way to express this in Sage? (I know this might be difficult as
python is dynamically typed).

This might be important for instance, when applying some operator to
the function. For instance: the results of one operation can be
different if the function is considered as one of a real variable, or
as one of a complex variable (say).

I think it would be important to have some unified framework for
functions, so that users
don't see that one function acts in one way and some other in a different way.

regards
Pablo

[Jason Grout]

Pablo De Napoli wrote:

> 2) the current implementation treats f(x)=sin(x) as a symbolic expression
>
> sage: f(x)=sin(x)
> sage: f
> x |--> sin(x)
> sage: type(f)
> <type 'sage.symbolic.expression.Expression'>
>
> However, mathematically a function is something different: for instance it
> has a domain and a range (say: this function takes one real number as
> argument, and returns another real number).
>
> �is there some way to express this in Sage? (I know this might be difficult as
> python is dynamically typed).
>
> This might be important for instance, when applying some operator to
> the function. For instance: the results of one operation can be
> different if the function is considered as one of a real variable, or
> as one of a complex variable (say).
>
> I think it would be important to have some unified framework for
> functions, so that users
> don't see that one function acts in one way and some other in a different way.

At one point, there was a discussion about being able to state the
domain and range of a function. Nicolas had some ideas about defining
that, but I can't find the post in which he talked about his ideas
(something about constructing domains and ranges like, say, the upper
half plane as Domain(RR^2,[y>0]) ). I'd love to be able to define a
function from R^n to R^m, for example, and have true vector-valued
functions.

I believe in 4.3 you can check the parents of something that is passed
into the evaluation function, so you can do some checking of the domain
and make the function respond differently depending on what type the
argument is.

[Robert Bradshaw]

On Dec 9, 2009, at 5:40 AM, Pablo De Napoli wrote:

> Many thanks to everybody for your help.
>
> some questions/remarks:
>
> 1) ¿Does every function needs to have two versions: a symbolic one and
> a numerical one?

Somewhat. The situation arises because one often only defines the
function numerically (i.e. it's an algorithm that takes numerical
input and returns numerical output), but wants to use it in symbolic
contexts.

> 2) the current implementation treats f(x)=sin(x) as a symbolic
> expression
>
> sage: f(x)=sin(x)
> sage: f
> x |--> sin(x)
> sage: type(f)
> <type 'sage.symbolic.expression.Expression'>
>
> However, mathematically a function is something different: for
> instance it
> has a domain and a range (say: this function takes one real number as
> argument, and returns another real number).
>
> ¿is there some way to express this in Sage? (I know this might be
> difficult as
> python is dynamically typed).
>
> This might be important for instance, when applying some operator to
> the function. For instance: the results of one operation can be
> different if the function is considered as one of a real variable, or
> as one of a complex variable (say).
>
> I think it would be important to have some unified framework for
> functions, so that users
> don't see that one function acts in one way and some other in a
> different way.

Right now functions do act differently depending on their input:

sage: log(2.0, 5.0)
0.430676558073393

sage: log(mod(2, 17), mod(5, 17))
6

If one looks at the source of log(a,b), the first thing it tries to do
is call a.log(b), which may do entirely different things depending on
what a is. I see this as a good thing--otherwise one would have to
have a whole bunch of log functions (log_real, log_complex,
descrete_log, ... (Actually, this last one would be an infinite
collection of log functions with distinct domains.)) That would
clutter the namespace and also make choosing the right one hard to
remember. Rather than view log as a function, I view it as a symbol
whose meaning is context sensitive.

Of course, there are situations where one would want to specify the
domain and range of a function, especially for user-defined functions.

- Robert