NaN (Not a Number) occurs in calculation of complex number for Bessel functions
DEAR PYTHON USERS DO MATHEMATICAL FUNCTIONS HAVE LIMITATION IN PYTHON in comparison with other programming languages ???? I have two mathematical functions: from scipy.special import sph_jn, sph_jnyn 1) sph_jn (n, z) ---> n is float, z is complex number for example: a,b=sph_jn ( 2.0 , 5+0.4j ) gives the following result:
a array( [ - 0.20416243 + 0.03963597j, - 0.10714653 - 0.06227716j, 0.13731305 - 0.07165432j ] ) b array( [ 0.10714653 + 0.06227716j, -0.15959617 + 0.06098154j, -0.18559289 - 0.01300886j ] )
2) sph_jnyn(n , x) --> n-float, x - float ,for example: c,d,e,f=sph_jnyn(2.0 , 3.0)
c array( [ 0.04704 , 0.3456775, 0.2986375 ] ) d array( [-0.3456775 , -0.18341166, 0.04704 ] ) e array( [ 0.3299975 , 0.06295916, -0.26703834 ] ) f array( [ -0.06295916, 0.28802472, 0.3299975 ] )
PROBLEM IS HERE!!! BUT , IF I GIVE ( it is necessary value for my program ): a , b = sph_jn ( 536 , 2513.2741228718346 + 201.0619298974676j ) I would like to see even very very deep comments as specific as possible!!!!!!
Happyman <bahtiyor_zohidov <at> mail.ru> writes: [clip]
IF I GIVE ( it is necessary value for my program ): a , b = sph_jn ( 536 , 2513.2741228718346 + 201.0619298974676j )
The implementation of the spherical Bessel functions is through this Fortran code: https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f#L... It does not have asymptotic expansions for dealing with parts of the complex plane where the computation via the recurrence does not work. -- Pauli Virtanen
Thanks Pauli But I have already very shortly built for bessel function, but the code you gave me is in Fortran.. I also used f2py but I could not manage to read fortran codes..that is why I have asked in Python what is wrong?? Пятница, 21 декабря 2012, 12:46 UTC от Pauli Virtanen <pav@iki.fi>:
Happyman <bahtiyor_zohidov <at> mail.ru> writes: [clip]
IF I GIVE ( it is necessary value for my program ): a , b = sph_jn ( 536 , 2513.2741228718346 + 201.0619298974676j )
The implementation of the spherical Bessel functions is through this Fortran code:
https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f#L...
It does not have asymptotic expansions for dealing with parts of the complex plane where the computation via the recurrence does not work.
-- Pauli Virtanen
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Happyman <bahtiyor_zohidov <at> mail.ru> writes:
Thanks Pauli But I have already very shortly built for bessel function, but the code you gave me is in Fortran.. I also used f2py but I could not manage to read fortran codes..that is why I have asked in Python what is wrong??
That Fortran code is `sph_jn`, which you used. It works using f2py. Only some of the special functions in scipy.special are written using Python as the language. Most of them are in C or in Fortran, using some existing special function library not written by us. Some of the implementations provided by these libraries are not complete, and do not cover the whole complex plane (or the real axis). Other functions (the more common ones), however, have very good implementations. -- Pauli Virtanen
I have everything in C or Fortran...According to my friends recommendations I started learning Python for my research... Do you mean the functions which gave Nan result has not been developed properly yet in Python, Don't you???? For about 1.5 months I have been facing the same problem for Bessel functions.. I think the code that I showed like an example is not working in Python. What to do ??? Пятница, 21 декабря 2012, 13:17 от Pauli Virtanen <pav@iki.fi>:
Happyman <bahtiyor_zohidov <at> mail.ru> writes:
Thanks Pauli But I have already very shortly built for bessel function, but the code you gave me is in Fortran.. I also used f2py but I could not manage to read fortran codes..that is why I have asked in Python what is wrong??
That Fortran code is `sph_jn`, which you used. It works using f2py.
Only some of the special functions in scipy.special are written using Python as the language. Most of them are in C or in Fortran, using some existing special function library not written by us.
Some of the implementations provided by these libraries are not complete, and do not cover the whole complex plane (or the real axis). Other functions (the more common ones), however, have very good implementations.
-- Pauli Virtanen
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On 12/21/2012 02:30 PM, Happyman wrote:
I have everything in C or Fortran...According to my friends recommendations I started learning Python for my research...
Do you mean the functions which gave Nan result has not been developed properly yet in Python, Don't you????
The way most of NumPy and SciPy works is by calling into C and Fortran code.
For about 1.5 months I have been facing the same problem for Bessel functions.. I think the code that I showed like an example is not working in Python. What to do ???
Do you have an implemention of the Bessel functions that work as you wish in C or Fortran? If so, that could be wrapped and called from Python. Dag Sverre
Dag Sverre Seljebotn <d.s.seljebotn <at> astro.uio.no> writes: [clip]
Do you have an implemention of the Bessel functions that work as you wish in C or Fortran? If so, that could be wrapped and called from Python.
For spherical Bessel functions it's possible to also use the relation to Bessel functions, which have a better implementation (AMOS) in Scipy: import numpy as np from scipy.special import jv def sph_jn(n, z): return jv(n + 0.5, z) * np.sqrt(np.pi / (2 * z)) print sph_jn(536, 2513.2741228718346 + 201.0619298974676j) # (2.5167666386507171e+81+3.3576357192536334e+81j) This should solve the problem. -- Pauli Virtanen
Received from Pauli Virtanen on Fri, Dec 21, 2012 at 08:59:02AM EST:
Dag Sverre Seljebotn <d.s.seljebotn <at> astro.uio.no> writes: [clip]
Do you have an implemention of the Bessel functions that work as you wish in C or Fortran? If so, that could be wrapped and called from Python.
For spherical Bessel functions it's possible to also use the relation to Bessel functions, which have a better implementation (AMOS) in Scipy:
import numpy as np from scipy.special import jv
def sph_jn(n, z): return jv(n + 0.5, z) * np.sqrt(np.pi / (2 * z))
print sph_jn(536, 2513.2741228718346 + 201.0619298974676j) # (2.5167666386507171e+81+3.3576357192536334e+81j)
This should solve the problem.
You can also try the spherical Bessel function in Sympy; it should be able to handle this case as well [1]. L.G. [1] http://docs.sympy.org/0.7.2/modules/functions/special.html#bessel-type-funct...
I think you advised about the code which is the same appearance. ========================================================================== Problem is not here Sir.... I will give you exactly what I was talking about. I have ready codes already(It would be kind of you if you checked the following codes, may be): ------------------------------------------------------------------------------------------------------------------------------ ## Bessel function of the first kind # mathematical form: Jn(x)--> n=arg1, x=arg2 # returns --> Jn(x) value in a complex form ------------------------------ def bes_1(arg1,arg2): nu=arg1+0.5 # nu=n+0.5: Jn(x) --> Jnu(x) return sqrt(pi/(2*arg2))*np.round(jv(nu,arg2),5) # jv(nu,arg2)--> from 'numpy.special' in PYTHON -------------------------------------------------------------------------------------------------------------------------------- ## Bessel function of the second kind # mathematical form: Yn(x)--> n=arg1, x=arg2 # returns --> Yn(x) value in a complex form --------------------------------- def bes_2 ( arg1, arg2 ): nu = arg1 + 0.5 # nu=n+0.5: Yn(x)--> Ynu(x) return sqrt ( pi / ( 2 * arg2 ) ) * np.round ( yv ( nu , arg2 ) , 5) # yv(nu,arg2)--> from 'numpy.special' in PYTHON --------------------------------- ------------------------------------------------------------------------------------------------------- ## Hankel function of the first kind # mathematical form: Hn(x)= Jn(x)+Yn(x)j # returns --> Hn(x) value in a complex form --------------------------------- def hank_1 ( arg1, arg2 ): return bes_1 ( arg1 , arg2 ) + bes_2 ( arg1 , arg2 ) * 1.0j # Hn(x)= Jn(x)+Yn(x)j ------------------------------------------------------------------------------------------------------------------------------------------- ## Bessel function of the first kind derivative # mathematical form: d(z*jn(z))=z*jn-1(z)-n*jn(z) where, z=x or m*x --------------------------------- def bes_der1 ( arg1 , arg2 ): return arg2 * bes_1 ( arg1 - 1, arg2 ) - arg1 * bes_1 ( arg1 , arg2 ) --------------------------------------------------------------------------------------------------------------------------------------------- ## Hankel function of the first kind derivative # mathematical form: d(z*hankeln(z))=z*hankeln-1(z)-n*hankeln(z) where, z=x or m*x def hank_der1 ( arg1 , arg2 ): return arg2 * hank_1 ( arg1 - 1, arg2 ) - arg1 * hank_1( arg1, arg2 ) --------------------------------------------------------------------------------------------------------------------------------------------- FOR MY CASE: m = 2513.2741228718346 + 201.0619298974676j x = 502.6548245743669 def F(m,x): nmax = x + 2.0 + 4.0 * x ** ( 1.0 / 3.0 ) # nmax= gives 536.0 as expected value nstop = np.round( nmax ) n = np.arange ( 0.0 ,nstop, dtype = float) z = m * x m2 = m * m val1 = m2 * bes_1 ( en , z ) * bes_der1 ( en, x) val2 = bes_1 ( en , x ) * bes_der1 ( en , z ) val3 = m2 * bes_1 ( en , z ) * hank_der1 ( en , x ) val4 = hank_1 ( en , x ) * bes_der1 ( en , z ) an = ( val1 - val2 ) / ( val3 - val4 ) val5 = bes_1 ( en , z ) * bes_der1 ( en, x ) val6 = bes_1 ( en , x ) * bes_der1 ( en, z ) val7 = bes_1 ( en , z ) * hank_der1 ( en, x ) val8 = hank_1 ( en , x ) * bes_der1 ( en , z ) bn = ( val5 - val6 ) / ( val7 - val8 ) return an, bn !!!! PROBLEM IS RETURNING THE an, bn at a given value which I showed before F(m,x) function WHAT IS WRONG WITH THIS?????? Пятница, 21 декабря 2012, 13:59 от Pauli Virtanen <pav@iki.fi>:
Dag Sverre Seljebotn <d.s.seljebotn <at> astro.uio.no> writes: [clip]
Do you have an implemention of the Bessel functions that work as you wish in C or Fortran? If so, that could be wrapped and called from Python.
For spherical Bessel functions it's possible to also use the relation to Bessel functions, which have a better implementation (AMOS) in Scipy:
import numpy as np from scipy.special import jv
def sph_jn(n, z): return jv(n + 0.5, z) * np.sqrt(np.pi / (2 * z))
print sph_jn(536, 2513.2741228718346 + 201.0619298974676j) # (2.5167666386507171e+81+3.3576357192536334e+81j)
This should solve the problem.
-- Pauli Virtanen
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi, Your code tries to to evaluate z = 1263309.3633394379 + 101064.74910119522j jv(536, z) # -> (inf+inf*j) In reality, this number is not infinite, but jv(536, z) == -2.3955170861527422e+43888 + 9.6910119847300024e+43887 These numbers (~ 10^43888) are too large for the floating point numbers that computers use (maximum ~ 10^308). This is why you get infinities and NaNs in the result. The same is true for the spherical Bessel functions. You will not be able to do this calculation using any software that uses only floating point numbers (Scipy, Matlab, ...). You need to use analytical properties of your problem to get rid of such large numbers. Alternatively, you can use arbitrary precision numbers. Python has libraries for that: http://code.google.com/p/mpmath/ By the way, the proper place for this discussion is the following mailing list: http://mail.scipy.org/mailman/listinfo/scipy-user -- Pauli Virtanen
Thanks But I could find for Win64 bit windows???? Second question: Did you mean that I have to put lens limits of those number??? Пятница, 21 декабря 2012, 15:45 UTC от Pauli Virtanen <pav@iki.fi>:
Hi,
Your code tries to to evaluate
z = 1263309.3633394379 + 101064.74910119522j jv(536, z) # -> (inf+inf*j)
In reality, this number is not infinite, but
jv(536, z) == -2.3955170861527422e+43888 + 9.6910119847300024e+43887
These numbers (~ 10^43888) are too large for the floating point numbers that computers use (maximum ~ 10^308). This is why you get infinities and NaNs in the result. The same is true for the spherical Bessel functions.
You will not be able to do this calculation using any software that uses only floating point numbers (Scipy, Matlab, ...).
You need to use analytical properties of your problem to get rid of such large numbers. Alternatively, you can use arbitrary precision numbers. Python has libraries for that:
http://code.google.com/p/mpmath/
By the way, the proper place for this discussion is the following mailing list:
http://mail.scipy.org/mailman/listinfo/scipy-user
-- Pauli Virtanen
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (4)
-
Dag Sverre Seljebotn -
Happyman -
Lev Givon -
Pauli Virtanen