# 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:

sage: var('r')
sage: bessel_J(0,r)


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.

edit retag close merge delete

Sort by » oldest newest most voted

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.

Here's a simple definition of a symbolic Bessel J:

----------------------------------------------------------------------
| 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:


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'

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

more

+1 for coming back and showing how it's done! I wish there were -- or wish I knew about! -- an equivalent Expression_factory.

( 2011-01-27 00:52:29 +0200 )edit

(1) Yep, it should return something symbolic, but unfortunately it doesn't right now. Some relevant tickets are #4102 and #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 ). It should be easy to make a symbolic wrapper around that. I've made this #10636.

more

OK, 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.

( 2011-01-16 20:18:26 +0200 )edit

FWIW, 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.

( 2011-01-16 21:17:42 +0200 )edit