# Ticket #3426: trac3426-fix-bessel-fns.patch

File trac3426-fix-bessel-fns.patch, 6.6 KB (added by AlexGhitza, 13 years ago)

apply instead of the previous patch

• ## sage/functions/special.py

# HG changeset patch
# User Alexandru Ghitza <aghitza@alum.mit.edu>
# Date 1222843216 -36000
# Node ID ed46740a68d3602b822352abf6378c70281b2b50
# Parent  4805d6f2365704ea92de52d00a97375299950add
trac 3426: fix Bessel functions in special.py

diff -r 4805d6f23657 -r ed46740a68d3 sage/functions/special.py
 a index (or "order") nu and argument z. INPUT: nu -- a real (or complex, for pari) number z  -- a real (positive) algorithm - "pari" or "maxima" or "scipy" prec - real precision (for Pari only) nu -- a real or complex number z  -- a real or complex number algorithm -- "pari" (default) or "maxima" or "scipy"; note that "maxima" and "scipy" can only deal with real nu and positive real z prec -- real precision (for Pari only) DEFINITION: \begin{verbatim} sage: bessel_I(1,1,"scipy") 0.565159103992... A corner case: sage: bessel_I(0,0) 1.00000000000000 If both nu and z are real, then bessel_I is known to be real if z is positive or if nu is an integer: sage: bessel_I(-2.1, 5.3) 22.5840990831073 sage: bessel_I(-2.1, -5.3) 21.4787545976446 - 6.97887041932783*I sage: bessel_I(-2, -5.3) 23.5424857046048 sage: bessel_I(10*I, 10) 1.19390297337705e6 - 688460.429969423*I """ if algorithm=="pari": from sage.libs.pari.all import pari R = RealField(prec) # this is a workaround for a bug in Pari 2.3.3 if nu.is_zero() and z.is_zero(): return R(1) C = ComplexField(prec) try: R = RealField(prec) nu = ZZ(nu) # pari doesn't like negative integers here, # but luckily I(-nu, z) = I(nu, z) if nu is integral nu = nu.abs() K = R except TypeError: K = C try: nu = R(nu) z = R(z) if (z > 0): K = R except TypeError: C = ComplexField(prec) nu = C(nu) z = C(z) K = C K = z.parent() return K(pari(nu).besseli(z, precision=prec)) pari_bes = pari(nu).besseli(z, precision=prec) if K is R: return K(pari_bes.real()) else: return K(pari_bes) elif algorithm=="scipy": if prec != 53: raise ValueError, "for the scipy algorithm the precision must be 53" 0.0583793793051868 sage: bessel_J(3,10,"scipy") 0.0583793793052... A corner case: sage: bessel_I(0,0) 1.00000000000000 If both nu and z are real, then bessel_J is known to be real if z is positive or if nu is an integer: sage: bessel_J(-2.1, 5.3) -0.123120872446240 sage: bessel_J(-2.1, -5.3) -0.117094908031941 + 0.0380464419481583*I sage: bessel_J(-2, -5.3) -0.0547481464525993 sage: bessel_J(10*I, 10) -119728.409042333 - 693848.987632644*I """ if algorithm=="pari": from sage.libs.pari.all import pari R = RealField(prec) # this is a workaround for a bug in Pari 2.3.3 if nu.is_zero() and z.is_zero(): return R(1) C = ComplexField(prec) fudge = 1 try: R = RealField(prec) nu = ZZ(nu) # pari doesn't like negative integers here, # but luckily J(-nu, z) = (-1)**nu J(nu, z) if nu is integral if nu < 0: nu = -nu fudge = (-1)**nu K = R except TypeError: K = C try: nu = R(nu) z = R(z) if (z > 0): K = R except TypeError: C = ComplexField(prec) nu = C(nu) z = C(z) K = C K = z.parent() return K(pari(nu).besselj(z, precision=prec)) pari_bes = pari(nu).besselj(z, precision=prec) if K is R: return fudge*K(pari_bes.real()) else: return fudge*K(pari_bes) elif algorithm=="scipy": if prec != 53: raise ValueError, "for the scipy algorithm the precision must be 53" otherwise. Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In Pari, nu can be complex and z must be real and positive. literature. In Pari, both nu and z can be complex; in SciPy they must be real and z mush be positive. EXAMPLES: sage: bessel_K(3,2,"scipy") 0.60 sage: bessel_K(1,1,"pari",100) 0.60190723019723457473754000154 sage: bessel_K(0,0) +Infinity If both nu and z are real, then bessel_K is known to be real if z is positive: sage: bessel_K(-2.1, 5.3) 0.00388720082599840 sage: bessel_K(-2.1, -5.3) 0.00369694767571369 - 70.9488385563182*I sage: bessel_K(-2, -5.3) 0.00375340283798982 - 73.9609001368291*I sage: bessel_K(10*I, 10) 9.82415743819925e-8 sage: bessel_K(2, -1.121) 1.23414145962983 - 0.547231658266306*I """ if algorithm=="scipy": if prec != 53: return sage_eval(ans) elif algorithm == 'pari': from sage.libs.pari.all import pari R = RealField(prec) C = ComplexField(prec) if z.is_zero() and nu.is_zero(): from sage.rings.infinity import Infinity return Infinity #TODO: sort out cases z = 0, nu != 0 try: R = RealField(prec) nu = R(nu) # K(-nu, z) = K(nu, z) for all nu nu = nu.abs() z = R(z) if (z > 0): K = R else: K = C except TypeError: C = ComplexField(prec) nu = C(nu) z = C(z) K = C K = z.parent() return K(pari(nu).besselk(z, precision=prec)) pari_bes = pari(nu).besselk(z, precision=prec) if K is R: return K(pari_bes.real()) else: return K(pari_bes) else: raise ValueError, "unknown algorithm '%s'"%algorithm