Ask Your Question
1

Bessel functions

asked 2011-01-16 09:44:24 -0500

Mike Witt gravatar image

updated 2011-06-16 08:46:12 -0500

Kelvin Li gravatar image

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 flag offensive close merge delete

2 answers

Sort by » oldest newest most voted
2

answered 2011-01-16 12:36:34 -0500

Mike Hansen gravatar image

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

edit flag offensive delete link more

Comments

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.

Mike Witt gravatar imageMike Witt ( 2011-01-16 13:18:26 -0500 )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.

Mike Witt gravatar imageMike Witt ( 2011-01-16 14:17:42 -0500 )edit
2

answered 2011-01-26 17:19:13 -0500

Mike Witt gravatar image

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()
edit flag offensive delete link more

Comments

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

DSM gravatar imageDSM ( 2011-01-26 17:52:29 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-01-16 09:44:24 -0500

Seen: 622 times

Last updated: Jan 26 '11