Ticket #15024: trac15024.patch

sage/functions/all.py

@@ -26 +26 @@
 
 from transcendental import (zeta, zetaderiv, zeta_symmetric, dickman_rho)
 
-from sage.functions.bessel import (bessel_I, bessel_J, bessel_K, bessel_Y, Bessel)
+from bessel import (Bessel, bessel_I, bessel_J, bessel_K, bessel_Y, hankel1,
+                    hankel2, spherical_bessel_J, spherical_bessel_Y,
+                    spherical_hankel1, spherical_hankel2)
 
 from special import (hypergeometric_U,
-                     spherical_bessel_J, spherical_bessel_Y,
-                     spherical_hankel1, spherical_hankel2,
                      spherical_harmonic, jacobi,
                      inverse_jacobi,
                      lngamma, error_fcn, elliptic_e,

sage/functions/bessel.py

@@ -90 +90 @@
    kind; they are also two linearly independent solutions of Bessel's
    differential equation. They are named for Hermann Hankel.
 
+-  When solving for separable solutions of Laplace's equation in
+   spherical coordinates, the radial equation has the form:
+
+   .. math::
+
+         x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0.
+
+
+   The spherical Bessel functions `j_n` and `y_n`,
+   are two linearly independent solutions to this equation. They are
+   related to the ordinary Bessel functions `J_n` and
+   `Y_n` by:
+
+   .. math::
+
+         j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x),
+
+   .. math::
+
+         y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x) = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x).
+
 EXAMPLES:
 
     Evaluate the Bessel J function symbolically and numerically::
@@ -143 +164 @@
 
 AUTHORS:
 
+    - Some of the documentation here has been adapted from David Joyner's
+      original documentation of Sage's special functions module (2006).
+
     - Benjamin Jones (2012-12-27): initial version
 
-    - Some of the documentation here has been adapted from David Joyner's
-      original documentation of Sage's special functions module (2006).
+    - Eviatar Bach (2013): Hankel and spherical Bessel and Hankel functions
 
 REFERENCES:
 
@@ -370 +393 @@
         EXAMPLES::
 
             sage: latex(bessel_J(1, x))
-            \operatorname{J_{1}}(x)
+            J_{1}(x)
         """
-        return r"\operatorname{J_{%s}}(%s)" % (latex(n), latex(z))
+        return r"J_{%s}(%s)" % (latex(n), latex(z))
 
 bessel_J = Function_Bessel_J()
 
@@ -530 +553 @@
         EXAMPLES::
 
             sage: latex(bessel_Y(1, x))
-            \operatorname{Y_{1}}(x)
+            Y_{1}(x)
         """
-        return r"\operatorname{Y_{%s}}(%s)" % (latex(n), latex(z))
+        return r"Y_{%s}(%s)" % (latex(n), latex(z))
 
 bessel_Y = Function_Bessel_Y()
 
@@ -703 +726 @@
         EXAMPLES::
 
             sage: latex(bessel_I(1, x))
-            \operatorname{I_{1}}(x)
+            I_{1}(x)
         """
-        return r"\operatorname{I_{%s}}(%s)" % (latex(n), latex(z))
+        return r"I_{%s}(%s)" % (latex(n), latex(z))
 
 bessel_I = Function_Bessel_I()
 
@@ -888 +911 @@
         EXAMPLES::
 
             sage: latex(bessel_K(1, x))
-            \operatorname{K_{1}}(x)
+            K_{1}(x)
         """
-        return r"\operatorname{K_{%s}}(%s)" % (latex(n), latex(z))
+        return r"K_{%s}(%s)" % (latex(n), latex(z))
 
 bessel_K = Function_Bessel_K()
 
@@ -1066 +1089 @@
     else:
         return _f
 
+
+class Hankel1(BuiltinFunction):
+    r"""
+    The Hankel function of the first kind
+
+    DEFINITION:
+
+    .. math::
+
+        H_\nu^{(1)}(z) = J_{\nu}(z) + iY_{\nu}(z)
+
+    EXAMPLES::
+
+        sage: hankel1(3, 4.)
+        0.430171473875622 - 0.182022115953485*I
+        sage: latex(hankel1(3, x))
+        H_{3}^{(1)}(x)
+        sage: hankel1(3., x).series(x == 2, 10).subs(x=3).n()  # abs tol 1e-12
+        0.309062682819597 - 0.512591541605233*I
+        sage: hankel1(3, 3.)
+        0.309062722255252 - 0.538541616105032*I
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'hankel1', nargs=2,
+                                 conversions=dict(maple='HankelH1',
+                                                  mathematica='HankelH1',
+                                                  maxima='hankel1',
+                                                  sympy='hankel1'))
+
+    def _eval_(self, nu, z):
+        nu, z = get_coercion_model().canonical_coercion(nu, z)
+        if is_inexact(nu) and not isinstance(nu, Expression):
+            return self._evalf_(nu, z, parent(nu))
+        return
+
+    def _evalf_(self, nu, z, parent):
+        from mpmath import hankel1
+        return mpmath_utils.call(hankel1, nu, z, parent=parent)
+
+    def _print_latex_(self, nu, z):
+        return r"H_{{{}}}^{{(1)}}({})".format(latex(nu), latex(z))
+
+    def _derivative_(self, nu, z, diff_param):
+        if diff_param == 1:
+            return (nu * hankel1(nu, z)) / z - hankel1(nu + 1, z)
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+hankel1 = Hankel1()
+
+
+class Hankel2(BuiltinFunction):
+    r"""
+    The Hankel function of the second kind
+
+    DEFINITION:
+
+    .. math::
+
+        H_\nu^{(2)}(z) = J_{\nu}(z) - iY_{\nu}(z)
+
+    EXAMPLES::
+
+        sage: hankel2(3, 4.)
+        0.430171473875622 + 0.182022115953485*I
+        sage: latex(hankel2(3, x))
+        H_{3}^{(2)}(x)
+        sage: hankel2(3., x).series(x == 2, 10).subs(x=3).n()  # abs tol 1e-12
+        0.309062682819597 + 0.512591541605234*I
+        sage: hankel2(3, 3.)
+        0.309062722255252 + 0.538541616105032*I
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'hankel2', nargs=2,
+                                 conversions=dict(maple='HankelH2',
+                                                  mathematica='HankelH2',
+                                                  maxima='hankel2',
+                                                  sympy='hankel2'))
+
+    def _eval_(self, nu, z):
+        nu, z = get_coercion_model().canonical_coercion(nu, z)
+        if is_inexact(nu) and not isinstance(nu, Expression):
+            return self._evalf_(nu, z, parent(nu))
+        return
+
+    def _evalf_(self, nu, z, parent):
+        from mpmath import hankel2
+        return mpmath_utils.call(hankel2, nu, z, parent=parent)
+
+    def _print_latex_(self, nu, z):
+        return r"H_{{{}}}^{{(2)}}({})".format(latex(nu), latex(z))
+
+    def _derivative_(self, nu, z, diff_param):
+        if diff_param == 1:
+            return (Integer(1) / 2) * (hankel2(nu - 1, z) - hankel2(nu + 1, z))
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+hankel2 = Hankel2()
+
+
+class SphericalBesselJ(BuiltinFunction):
+    r"""
+    The spherical Bessel function of the first kind
+
+    DEFINITION:
+
+    .. math::
+
+        j_n(z) = \sqrt{\frac{1}{2}\pi/z} J_{n + \frac{1}{2}}(z)
+
+    EXAMPLES::
+
+        sage: spherical_bessel_J(3., x).series(x == 2, 10).subs(x=3).n()
+        0.152051648665037
+        sage: spherical_bessel_J(3., 3)
+        0.152051662030533
+        sage: spherical_bessel_J(4, x).simplify()
+        -((45/x^2 - 105/x^4 - 1)*sin(x) + 5*(21/x^2 - 2)*cos(x)/x)/x
+        sage: latex(spherical_bessel_J(4, x))
+        j_{4}(x)
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'spherical_bessel_J', nargs=2,
+                                 conversions=dict(mathematica=
+                                                  'SphericalBesselJ',
+                                                  maxima='spherical_bessel_j',
+                                                  sympy='jn'))
+
+    def _eval_(self, n, z):
+        n, z = get_coercion_model().canonical_coercion(n, z)
+        if is_inexact(n) and not isinstance(n, Expression):
+            return self._evalf_(n, z, parent(n))
+        return
+
+    def _evalf_(self, n, z, parent):
+        return mpmath_utils.call(spherical_bessel_f, 'besselj', n, z,
+                                 parent=parent)
+
+    def _print_latex_(self, n, z):
+        return r"j_{{{}}}({})".format(latex(n), latex(z))
+
+    def _derivative_(self, n, z, diff_param):
+        if diff_param == 1:
+            return (spherical_bessel_J(n - 1, z) -
+                    ((n + 1) / z) * spherical_bessel_J(n, z))
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+spherical_bessel_J = SphericalBesselJ()
+
+
+class SphericalBesselY(BuiltinFunction):
+    r"""
+    The spherical Bessel function of the second kind
+
+    DEFINITION:
+
+    .. math::
+
+        y_n(z) = \sqrt{\frac{1}{2}\pi/z} Y_{n + \frac{1}{2}}(z)
+
+    EXAMPLES::
+
+        sage: spherical_bessel_Y(-3, x).simplify()
+        ((3/x^2 - 1)*sin(x) - 3*cos(x)/x)/x
+        sage: spherical_bessel_Y(3 + 2 * I, 5 - 0.2 * I)
+        -0.270205813266440 - 0.615994702714957*I
+        sage: integrate(spherical_bessel_Y(0, x), x)
+        -1/2*Ei(I*x) - 1/2*Ei(-I*x)
+        sage: latex(spherical_bessel_Y(0, x))
+        y_{0}(x)
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'spherical_bessel_Y', nargs=2,
+                                 conversions=dict(mathematica=
+                                                  'SphericalBesselY',
+                                                  maxima='spherical_bessel_y',
+                                                  sympy='yn'))
+
+    def _eval_(self, n, z):
+        n, z = get_coercion_model().canonical_coercion(n, z)
+        if is_inexact(n) and not isinstance(n, Expression):
+            return self._evalf_(n, z, parent(n))
+        return
+
+    def _evalf_(self, n, z, parent):
+        return mpmath_utils.call(spherical_bessel_f, 'bessely', n, z,
+                                 parent=parent)
+
+    def _print_latex_(self, n, z):
+        return r"y_{{{}}}({})".format(latex(n), latex(z))
+
+    def _derivative_(self, n, z, diff_param):
+        if diff_param == 1:
+            return (-spherical_bessel_Y(n, z) / (2 * z) +
+                    (spherical_bessel_Y(n - 1, z) -
+                     spherical_bessel_Y(n + 1, z)) / 2)
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+spherical_bessel_Y = SphericalBesselY()
+
+
+class SphericalHankel1(BuiltinFunction):
+    r"""
+    The spherical Hankel function of the first kind
+
+    DEFINITION:
+
+    .. math::
+
+        h_n^{(1)}(z) = \sqrt{\frac{1}{2}\pi/z} H_{n + \frac{1}{2}}^{(1)}(z)
+
+    EXAMPLES::
+
+        sage: spherical_hankel1(1, x).simplify()
+        -(x + I)*e^(I*x)/x^2
+        sage: spherical_hankel1(3 + 2 * I, 5 - 0.2 * I)
+        1.25375216869913 - 0.518011435921789*I
+        sage: integrate(spherical_hankel1(3, x), x)
+        Ei(I*x) - 6*gamma(-1, -I*x) - 15*gamma(-2, -I*x) - 15*gamma(-3, -I*x)
+        sage: latex(spherical_hankel1(3, x))
+        h_{3}^{(1)}(x)
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'spherical_hankel1', nargs=2,
+                                 conversions=dict(mathematica=
+                                                  'SphericalHankelH1',
+                                                  maxima='spherical_hankel1'))
+
+    def _eval_(self, n, z):
+        n, z = get_coercion_model().canonical_coercion(n, z)
+        if is_inexact(n) and not isinstance(n, Expression):
+            return self._evalf_(n, z, parent(n))
+        return
+
+    def _evalf_(self, n, z, parent):
+        return mpmath_utils.call(spherical_bessel_f, 'hankel1', n, z,
+                                 parent=parent)
+
+    def _print_latex_(self, n, z):
+        return r"h_{{{}}}^{{(1)}}({})".format(latex(n), latex(z))
+
+    def _derivative_(self, n, z, diff_param):
+        if diff_param == 1:
+            return (-spherical_hankel1(n, z) / (2 * z) +
+                    (spherical_hankel1(n - 1, z) -
+                     spherical_hankel1(n + 1, z)) / 2)
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+spherical_hankel1 = SphericalHankel1()
+
+
+class SphericalHankel2(BuiltinFunction):
+    r"""
+    The spherical Hankel function of the second kind
+
+    DEFINITION:
+
+    .. math::
+
+        h_n^{(2)}(z) = \sqrt{\frac{1}{2}\pi/z} H_{n + \frac{1}{2}}^{(2)}(z)
+
+    EXAMPLES::
+
+        sage: spherical_hankel2(1, x).simplify()
+        -(x - I)*e^(-I*x)/x^2
+        sage: spherical_hankel2(3 + 2*I, 5 - 0.2*I)
+        0.0217627632692163 + 0.0224001906110906*I
+        sage: integrate(spherical_hankel2(3, x), x)
+        Ei(-I*x) - 6*gamma(-1, I*x) - 15*gamma(-2, I*x) - 15*gamma(-3, I*x)
+        sage: latex(spherical_hankel2(3, x))
+        h_{3}^{(2)}(x)
+    """
+    def __init__(self):
+        BuiltinFunction.__init__(self, 'spherical_hankel2', nargs=2,
+                                 conversions=dict(mathematica=
+                                                  'SphericalHankelH2',
+                                                  maxima='spherical_hankel2'))
+
+    def _eval_(self, n, z):
+        n, z = get_coercion_model().canonical_coercion(n, z)
+        if is_inexact(n) and not isinstance(n, Expression):
+            return self._evalf_(n, z, parent(n))
+        return
+
+    def _evalf_(self, n, z, parent):
+        return mpmath_utils.call(spherical_bessel_f, 'hankel2', n, z,
+                                 parent=parent)
+
+    def _print_latex_(self, n, z):
+        return r"h_{{{}}}^{{(2)}}({})".format(latex(n), latex(z))
+
+    def _derivative_(self, n, z, diff_param):
+        if diff_param == 1:
+            return (-spherical_hankel2(n, z) / (2 * z) +
+                    (spherical_hankel2(n - 1, z) -
+                     spherical_hankel2(n + 1, z)) / 2)
+        else:
+            raise NotImplementedError('derivative with respect to order')
+
+spherical_hankel2 = SphericalHankel2()
+
+
+def spherical_bessel_f(F, n, z):
+    r"""
+    Numerically evaluate the spherical version, `f`, of the Bessel function `F`
+    by computing `f_n(z) = \sqrt{\frac{1}{2}\pi/z} F_{n + \frac{1}{2}}(z)`.
+    According to Abramowitz & Stegun, this identity holds for the Bessel
+    functions `J`, `Y`, `K`, `I`, `H^{(1)}`, and `H^{(2)}`.
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import spherical_bessel_f
+        sage: spherical_bessel_f('besselj', 3, 4)
+        mpf('0.22924385795503024')
+        sage: spherical_bessel_f('hankel1', 3, 4)
+        mpc(real='0.22924385795503024', imag='-0.21864196590306359')
+    """
+    from mpmath import mp
+    ctx = mp
+    prec = ctx.prec
+    try:
+        n = ctx.convert(n)
+        z = ctx.convert(z)
+        ctx.prec += 10
+        Fz = getattr(ctx, F)(n + 0.5, z)
+        hpi = 0.5 * ctx.pi()
+        ctx.prec += 10
+        hpioz = hpi / z
+        ctx.prec += 10
+        sqrthpioz = ctx.sqrt(hpioz)
+        ctx.prec += 10
+        return sqrthpioz * Fz
+    finally:
+        ctx.prec = prec
+
 ####################################################
 ###  to be removed after the deprecation period  ###
 ####################################################

sage/functions/special.py

@@ -105 +105 @@
    
    .. math::
 
-     \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta . 
-
-
-
--  When solving for separable solutions of Laplace's equation in
-   spherical coordinates, the radial equation has the form:
-   
-   .. math::
-
-         x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0. 
-
-
-   The spherical Bessel functions `j_n` and `y_n`,
-   are two linearly independent solutions to this equation. They are
-   related to the ordinary Bessel functions `J_n` and
-   `Y_n` by:
-   
-   .. math::
-
-         j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x), 
-
-
-   
-   .. math::
-
-         y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x)     = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x).  
-
-
+     \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta .
 
 -  For `x>0`, the confluent hypergeometric function
    `y = U(a,b,x)` is defined to be the solution to Kummer's
@@ -666 +639 @@
     else:
         raise ValueError, "unknown algorithm '%s'"%algorithm
 
-def spherical_bessel_J(n, var, algorithm="maxima"):
-    r"""
-    Returns the spherical Bessel function of the first kind for
-    integers n >= 1.
-    
-    Reference: AS 10.1.8 page 437 and AS 10.1.15 page 439.
-    
-    EXAMPLES::
-    
-        sage: spherical_bessel_J(2,x)
-        ((3/x^2 - 1)*sin(x) - 3*cos(x)/x)/x
-        sage: spherical_bessel_J(1, 5.2, algorithm='scipy')
-        -0.12277149950007...
-        sage: spherical_bessel_J(1, 3, algorithm='scipy')
-        0.345677499762355...
-    """
-    if algorithm=="scipy":
-        from scipy.special.specfun import sphj
-        return sphj(int(n), float(var))[1][-1]
-    elif algorithm == 'maxima':
-        _init()
-        return meval("spherical_bessel_j(%s,%s)"%(ZZ(n),var))
-    else:
-        raise ValueError, "unknown algorithm '%s'"%algorithm
-    
-def spherical_bessel_Y(n,var, algorithm="maxima"):
-    r"""
-    Returns the spherical Bessel function of the second kind for
-    integers n -1.
-    
-    Reference: AS 10.1.9 page 437 and AS 10.1.15 page 439.
-    
-    EXAMPLES::
-    
-        sage: x = PolynomialRing(QQ, 'x').gen()
-        sage: spherical_bessel_Y(2,x)
-        -((3/x^2 - 1)*cos(x) + 3*sin(x)/x)/x
-    """
-    if algorithm=="scipy":
-        import scipy.special
-        ans = str(scipy.special.sph_yn(int(n),float(var)))
-        ans = ans.replace("(","")
-        ans = ans.replace(")","")
-        ans = ans.replace("j","*I")
-        return sage_eval(ans)
-    elif algorithm == 'maxima':
-        _init()
-        return meval("spherical_bessel_y(%s,%s)"%(ZZ(n),var))
-    else:
-        raise ValueError, "unknown algorithm '%s'"%algorithm
-
-def spherical_hankel1(n, var):
-    r"""
-    Returns the spherical Hankel function of the first kind for
-    integers `n > -1`, written as a string. Reference: AS
-    10.1.36 page 439.
-    
-    EXAMPLES::
-    
-        sage: spherical_hankel1(2, x)
-        (I*x^2 - 3*x - 3*I)*e^(I*x)/x^3
-    """
-    return maxima_function("spherical_hankel1")(ZZ(n), var)
-
-def spherical_hankel2(n,x):
-    r"""
-    Returns the spherical Hankel function of the second kind for
-    integers `n > -1`, written as a string. Reference: AS 10.1.17 page
-    439.
-    
-    EXAMPLES::
-    
-        sage: spherical_hankel2(2, x)
-        (-I*x^2 - 3*x + 3*I)*e^(-I*x)/x^3
-    
-    Here I = sqrt(-1).
-    """
-    return maxima_function("spherical_hankel2")(ZZ(n), x)
-
 def spherical_harmonic(m,n,x,y):
     r"""
     Returns the spherical Harmonic function of the second kind for