# Ticket #15024: trac15024.patch

File trac15024.patch, 19.1 KB (added by eviatarbach, 8 years ago)
• ## sage/functions/all.py

# HG changeset patch
# User Eviatar Bach <eviatarbach@gmail.com>
# Date 1376006015 25200
# Node ID 684d2ce97584d9c2424306687cd94b4c0937b4b6
diff --git a/sage/functions/all.py b/sage/functions/all.py
diff --git a/sage/functions/bessel.py b/sage/functions/bessel.py
 a 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:: 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: 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() 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() 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() 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() 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  ### ####################################################
diff --git a/sage/functions/special.py b/sage/functions/special.py
 a .. 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 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