Bessel functions

edit

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