Hi,
for testing that my FEM calculation of the radial Schroedinger
equation is correct, I test the claculated eigenvalues. One example is
a 3D spherical well, the energies of which are zeros of the spherical
bessel function. Here is how to calculate in sympy (resp. mpmath):
from sympy.mpmath import *
def jn(n, x):
return sqrt(pi/2*x) * besselj(n+0.5, x)
f = lambda x: jn(0, x)
nprint([findroot(f, k) for k in [3, 6, 9, 12]])
f = lambda x: jn(1, x)
nprint([findroot(f, k) for k in [4, 7, 10, 13]])
f = lambda x: jn(2, x)
nprint([findroot(f, k) for k in [5, 9, 12, 15]])
f = lambda x: jn(3, x)
nprint([findroot(f, k) for k in [6, 10, 13, 16]])
f = lambda x: jn(4, x)
nprint([findroot(f, k) for k in [8, 11, 15, 18]])
f = lambda x: jn(5, x)
nprint([findroot(f, k) for k in [9, 12, 16, 19]])
this prints:
[3.14159, 6.28319, 9.42478, 12.5664]
[4.49341, 7.72525, 10.9041, 14.0662]
[5.76346, 9.09501, 12.3229, 15.5146]
[6.98793, 10.4171, 13.698, 16.9236]
[8.18256, 11.7049, 15.0397, 18.3013]
[9.35581, 12.9665, 16.3547, 19.6532]
which is correct, as can be checked with the table at:
http://quantummechanics.ucsd.edu/ph130a/130_notes/node226.html
We should in my opinion have the spherical bessel functions in sympy.
There are bessel functions in mpmath and they are related to the
spherical ones easily, see the function jn(n, x) above. Mathematica
has them too:
http://reference.wolfram.com/mathematica/ref/SphericalBesselJ.html
We should have a function bessel_j_zero, like mathematica:
http://mathworld.wolfram.com/BesselFunctionZeros.html
As far as I know mathematica doesn't have spherical_bessel_j_zero, but
in my opinion it will be handy.
What is the best way to do that? It seems to me the distance between
the roots is approximately pi, so one algorithm that works for me is:
def jn_zero(n, k):
"""
Calculates the k-th zero of the spherical bessel function j_n.
"""
f = lambda x: jn(n, x)
# we need to approximate the position of the first root:
root = n+pi
# determine the first root exactly:
root = findroot(f, root)
for i in range(k-1):
# estimate the position of the next root using the last root + pi:
root = findroot(f, root+pi)
return root
nprint([jn_zero(0, k) for k in [1, 2, 3, 4]])
nprint([jn_zero(1, k) for k in [1, 2, 3, 4]])
nprint([jn_zero(2, k) for k in [1, 2, 3, 4]])
nprint([jn_zero(3, k) for k in [1, 2, 3, 4]])
nprint([jn_zero(4, k) for k in [1, 2, 3, 4]])
nprint([jn_zero(5, k) for k in [1, 2, 3, 4]])
This exactly reproduces the table above. It should probably cache the values.
Ondrej
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com
To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
