ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 26 Jan 2011 17:52:29 -0600Bessel functionshttp://ask.sagemath.org/question/7878/bessel-functions/A couple of questions about Bessel functions:
(1) Shouldn't this return return a symbolic expression, so that it can be used as part of an expression to be evaluated later with a specific r:<pre>
sage: var('r')
sage: bessel_J(0,r)
</pre>
I'm thinking specifically of building up a series that will later be plotted.
(2) There was talk a while back about making the zeros of the Bessel functions available without having to calculate them with find_root() (Not that there's a problem doing that, but it would just be a convenience.) I was just wondering if anything came of that.
Sun, 16 Jan 2011 09:44:24 -0600http://ask.sagemath.org/question/7878/bessel-functions/Answer by Mike Hansen for <p>A couple of questions about Bessel functions:</p>
<p>(1) Shouldn't this return return a symbolic expression, so that it can be used as part of an expression to be evaluated later with a specific r:</p><pre>sage: var('r')
sage: bessel_J(0,r)
</pre><p></p>
<p>I'm thinking specifically of building up a series that will later be plotted.</p>
<p>(2) There was talk a while back about making the zeros of the Bessel functions available without having to calculate them with find_root() (Not that there's a problem doing that, but it would just be a convenience.) I was just wondering if anything came of that.</p>
http://ask.sagemath.org/question/7878/bessel-functions/?answer=11979#post-id-11979(1) Yep, it should return something symbolic, but unfortunately it doesn't right now. Some relevant tickets are [#4102](http://trac.sagemath.org/sage_trac/ticket/4102) and [#7636](http://trac.sagemath.org/sage_trac/ticket/7636) (see the referenced thread). I'll try to make it a priority to get that done for 4.6.2.
(2) I don't think anything's been done with that yet either. mpmath provides efficient numerical evaluation of the zeros (see [here](http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html) ). It should be easy to make a symbolic wrapper around that. I've made this [#10636](http://trac.sagemath.org/sage_trac/ticket/10636).Sun, 16 Jan 2011 12:36:34 -0600http://ask.sagemath.org/question/7878/bessel-functions/?answer=11979#post-id-11979Comment by Mike Witt for <p>(1) Yep, it should return something symbolic, but unfortunately it doesn't right now. Some relevant tickets are <a href="http://trac.sagemath.org/sage_trac/ticket/4102">#4102</a> and <a href="http://trac.sagemath.org/sage_trac/ticket/7636">#7636</a> (see the referenced thread). I'll try to make it a priority to get that done for 4.6.2.</p>
<p>(2) I don't think anything's been done with that yet either. mpmath provides efficient numerical evaluation of the zeros (see <a href="http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html">here</a> ). It should be easy to make a symbolic wrapper around that. I've made this <a href="http://trac.sagemath.org/sage_trac/ticket/10636">#10636</a>.</p>
http://ask.sagemath.org/question/7878/bessel-functions/?comment=22273#post-id-22273OK, thanks. Like I mentioned #2 is just a convenience. But #1 really had me stumped. I do see that there appears to be an explanation of how to create a symbolic function in one of those threads, and I'll try that out.Sun, 16 Jan 2011 13:18:26 -0600http://ask.sagemath.org/question/7878/bessel-functions/?comment=22273#post-id-22273Comment by Mike Witt for <p>(1) Yep, it should return something symbolic, but unfortunately it doesn't right now. Some relevant tickets are <a href="http://trac.sagemath.org/sage_trac/ticket/4102">#4102</a> and <a href="http://trac.sagemath.org/sage_trac/ticket/7636">#7636</a> (see the referenced thread). I'll try to make it a priority to get that done for 4.6.2.</p>
<p>(2) I don't think anything's been done with that yet either. mpmath provides efficient numerical evaluation of the zeros (see <a href="http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html">here</a> ). It should be easy to make a symbolic wrapper around that. I've made this <a href="http://trac.sagemath.org/sage_trac/ticket/10636">#10636</a>.</p>
http://ask.sagemath.org/question/7878/bessel-functions/?comment=22272#post-id-22272FWIW, I can't seem to get the barebones example of a symbolic bessel_J from the above thread working. But I do look forward to any progress on those two tickets. They seem to apply far beyond this one particular issue.Sun, 16 Jan 2011 14:17:42 -0600http://ask.sagemath.org/question/7878/bessel-functions/?comment=22272#post-id-22272Answer by Mike Witt for <p>A couple of questions about Bessel functions:</p>
<p>(1) Shouldn't this return return a symbolic expression, so that it can be used as part of an expression to be evaluated later with a specific r:</p><pre>sage: var('r')
sage: bessel_J(0,r)
</pre><p></p>
<p>I'm thinking specifically of building up a series that will later be plotted.</p>
<p>(2) There was talk a while back about making the zeros of the Bessel functions available without having to calculate them with find_root() (Not that there's a problem doing that, but it would just be a convenience.) I was just wondering if anything came of that.</p>
http://ask.sagemath.org/question/7878/bessel-functions/?answer=12017#post-id-12017Just as a follow up to this. Since, even though creating a symbolic function is "easy", it wasn't easy (for me) to locate the information on how to do it. So just in case someone else might benefit from this, I thought I'd post the following examples which work on Sage 4.6.
Here's a simple definition of a symbolic Bessel J:
<pre>
----------------------------------------------------------------------
| Sage Version 4.6, Release Date: 2010-10-30 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: from sage.symbolic.function_factory import function_factory
sage: def j_eval(self, m, z, parent=None): return(bessel_J(m,z))
....:
sage: BesselJ = function_factory( name='BesselJ', nargs=2, evalf_func=j_eval)
sage:
sage: var('m,z')
(m, z)
sage: BesselJ(m,z)
BesselJ(m, z)
sage:
sage: integral_numerical(BesselJ(0,z),0,1)
(0.91973041008976042, 1.0211058777295894e-14)
sage:
</pre>
And here's the code I'm actually currently using. For whatever that's worth. Any comments on how to improve it would certainly be appreciated. The 'zero' table could obviously be made 'growable'
<pre>
from sage.symbolic.function_factory import function_factory
class Bessel:
def __init__(self):
self.JZ_max_m = -1
self.JZ_max_n = -1
self.TSIZE = 10
def sJ_evalf(self, m, z, parent=None): return(bessel_J(m,z))
def sJ_latex(self, m, z): return('\\mathbf{J}_{m}\\left(z\\right)')
self.J = function_factory(name='Bessel_J',
nargs=2, print_latex_func=sJ_latex, evalf_func=sJ_evalf)
#def J(self, m, n): return( bessel_J(m,n) )
def I(self, m, n): return( bessel_I(m,n) )
# For now either we have the table or we don't. But in the future I
# might like to be able to "extend" it, starting from what I've got.
def JZ(self, m, n):
if m > self.JZ_max_m or n > self.JZ_max_n:
print('Creating new table for Bessel.JZ ...')
self.JZ_table = []
for mm in range(self.TSIZE):
JZ_row = []
def jn(z): return bessel_J(mm, z)
low = 0; high = 1
for nn in range(self.TSIZE):
root = -1
while root == -1:
try: root = find_root(lambda x:jn(x), low, high)
except (RuntimeError):
low = high
high = low+1
continue
#print('Found root: %s' %root)
low = root+1; high = low+1
#print('New interval: %s to %s' %(low,high))
JZ_row += [root,]
self.JZ_table += [JZ_row,]
self.JZ_max_m = self.TSIZE
self.JZ_max_n = self.TSIZE
print('Done')
return(self.JZ_table[m][n])
def Bessel_Test():
p1 = plot(lambda x:bessel_J(0, x), (x, 0, 32), linestyle='--')
b = Bessel()
p1 += plot(lambda x:b.J(0, x), (x, 0, 32), color='yellow')
for r in range(10):
p1 += point([b.JZ(0,r), 0])
p1.show()
p2 = plot(lambda x:bessel_J(9, x), (x, 0, 42), linestyle='--')
p2 += plot(lambda x:b.J(9, x), (x, 0, 42), color='pink')
for r in range(10):
p2 += point([b.JZ(9,r), 0])
p2.show()
</pre>Wed, 26 Jan 2011 17:19:13 -0600http://ask.sagemath.org/question/7878/bessel-functions/?answer=12017#post-id-12017Comment by DSM for <p>Just as a follow up to this. Since, even though creating a symbolic function is "easy", it wasn't easy (for me) to locate the information on how to do it. So just in case someone else might benefit from this, I thought I'd post the following examples which work on Sage 4.6.</p>
<p>Here's a simple definition of a symbolic Bessel J:</p>
<pre>----------------------------------------------------------------------
| Sage Version 4.6, Release Date: 2010-10-30 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: from sage.symbolic.function_factory import function_factory
sage: def j_eval(self, m, z, parent=None): return(bessel_J(m,z))
....:
sage: BesselJ = function_factory( name='BesselJ', nargs=2, evalf_func=j_eval)
sage:
sage: var('m,z')
(m, z)
sage: BesselJ(m,z)
BesselJ(m, z)
sage:
sage: integral_numerical(BesselJ(0,z),0,1)
(0.91973041008976042, 1.0211058777295894e-14)
sage:
</pre>
<p>And here's the code I'm actually currently using. For whatever that's worth. Any comments on how to improve it would certainly be appreciated. The 'zero' table could obviously be made 'growable'</p>
<pre>from sage.symbolic.function_factory import function_factory
class Bessel:
def __init__(self):
self.JZ_max_m = -1
self.JZ_max_n = -1
self.TSIZE = 10
def sJ_evalf(self, m, z, parent=None): return(bessel_J(m,z))
def sJ_latex(self, m, z): return('\\mathbf{J}_{m}\\left(z\\right)')
self.J = function_factory(name='Bessel_J',
nargs=2, print_latex_func=sJ_latex, evalf_func=sJ_evalf)
#def J(self, m, n): return( bessel_J(m,n) )
def I(self, m, n): return( bessel_I(m,n) )
# For now either we have the table or we don't. But in the future I
# might like to be able to "extend" it, starting from what I've got.
def JZ(self, m, n):
if m > self.JZ_max_m or n > self.JZ_max_n:
print('Creating new table for Bessel.JZ ...')
self.JZ_table = []
for mm in range(self.TSIZE):
JZ_row = []
def jn(z): return bessel_J(mm, z)
low = 0; high = 1
for nn in range(self.TSIZE):
root = -1
while root == -1:
try: root = find_root(lambda x:jn(x), low, high)
except (RuntimeError):
low = high
high = low+1
continue
#print('Found root: %s' %root)
low = root+1; high = low+1
#print('New interval: %s to %s' %(low,high))
JZ_row += [root,]
self.JZ_table += [JZ_row,]
self.JZ_max_m = self.TSIZE
self.JZ_max_n = self.TSIZE
print('Done')
return(self.JZ_table[m][n])
def Bessel_Test():
p1 = plot(lambda x:bessel_J(0, x), (x, 0, 32), linestyle='--')
b = Bessel()
p1 += plot(lambda x:b.J(0, x), (x, 0, 32), color='yellow')
for r in range(10):
p1 += point([b.JZ(0,r), 0])
p1.show()
p2 = plot(lambda x:bessel_J(9, x), (x, 0, 42), linestyle='--')
p2 += plot(lambda x:b.J(9, x), (x, 0, 42), color='pink')
for r in range(10):
p2 += point([b.JZ(9,r), 0])
p2.show()
</pre>
http://ask.sagemath.org/question/7878/bessel-functions/?comment=22236#post-id-22236+1 for coming back and showing how it's done! I wish there were -- or wish I knew about! -- an equivalent Expression_factory.Wed, 26 Jan 2011 17:52:29 -0600http://ask.sagemath.org/question/7878/bessel-functions/?comment=22236#post-id-22236