Ticket #4102: trac_symbolic_bessel_v7.2.patch

File trac_symbolic_bessel_v7.2.patch, 74.2 KB (added by benjaminfjones, 8 years ago)
  • doc/en/reference/functions/index.rst

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1363120860 25200
    # Node ID 101556ef73af6f4eec47e5e9b9cba00e5ae1579d
    # Parent  865ef62dff12b716149e193c3cb62bc8af126ffc
    Trac 4102: Implement symbolic Bessel functions
    
    diff --git a/doc/en/reference/functions/index.rst b/doc/en/reference/functions/index.rst
    a b  
    1111   sage/functions/orthogonal_polys
    1212   sage/functions/other
    1313   sage/functions/special
     14   sage/functions/bessel
    1415   sage/functions/exp_integral
    1516   sage/functions/wigner
    1617   sage/functions/generalized
  • sage/functions/all.py

    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    2626
    2727from transcendental import (zeta, zetaderiv, zeta_symmetric, dickman_rho)
    2828
    29 from special import (bessel_I, bessel_J, bessel_K, bessel_Y,
    30                      hypergeometric_U, Bessel,
     29from sage.functions.bessel import (bessel_I, bessel_J, bessel_K, bessel_Y, Bessel)
     30
     31from special import (hypergeometric_U,
    3132                     spherical_bessel_J, spherical_bessel_Y,
    3233                     spherical_hankel1, spherical_hankel2,
    3334                     spherical_harmonic, jacobi,
     
    3637                     elliptic_f, elliptic_ec, elliptic_eu,
    3738                     elliptic_kc, elliptic_pi, elliptic_j,
    3839                     airy_ai, airy_bi)
    39                        
     40
    4041from orthogonal_polys import (chebyshev_T,
    4142                              chebyshev_U,
    4243                              gen_laguerre,
  • new file sage/functions/bessel.py

    diff --git a/sage/functions/bessel.py b/sage/functions/bessel.py
    new file mode 100644
    - +  
     1r"""
     2Bessel Functions
     3
     4This module provides symbolic Bessel Functions. These functions use the
     5`mpmath Library`_ for numerical evaluation and Maxima, GiNaC, Pynac for
     6symbolics.
     7
     8The main objects which are exported from this module are:
     9
     10 * `bessel_J` - The Bessel J function
     11 * `bessel_Y` - The Bessel J function
     12 * `bessel_I` - The Bessel J function
     13 * `bessel_K` - The Bessel J function
     14 * `Bessel`   - A factory function for producing Bessel functions of
     15   various kinds and orders.
     16
     17-  Bessel functions, first defined by the Swiss mathematician
     18   Daniel Bernoulli and named after Friedrich Bessel, are canonical
     19   solutions y(x) of Bessel's differential equation:
     20
     21   .. math::
     22
     23         x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \nu^2)y = 0,
     24
     25   for an arbitrary real number `\nu` (the order).
     26
     27-  In this module, `J_\nu` denotes the unique solution of Bessel's equation
     28   which is non-singular at `x = 0`. This function is known as the Bessel
     29   Function of the First Kind. This function also arises as a special case
     30   of the hypergeometric function `{}_0F_1`:
     31
     32   .. math::
     33
     34        J_\nu(x) = \frac{x^n}{2^\nu \Gamma(\nu + 1)} {}_0F_1(\nu +
     35        1, -\frac{x^2}{4}).
     36
     37-  The second linearly independent solution to Bessel's equation (which is
     38   singular at `x=0`) is denoted by `Y_\nu` and is called the Bessel
     39   Function of the Second Kind:
     40
     41   .. math::
     42
     43        Y_\nu(x) = \frac{ J_\nu(x) \cos(\pi \nu) -
     44        J_{-\nu}(x)}{\sin(\pi \nu)}.
     45
     46-  There are also two commonly used combinations of the Bessel J and Y
     47   Functions. The Bessel I Function, or the Modified Bessel Function of the
     48   First Kind, is defined by:
     49
     50   .. math::
     51
     52       I_\nu(x) = i^{-\nu} J_\nu(ix).
     53
     54   The Bessel K Function, or the Modified Bessel Function of the Second Kind,
     55   is defined by:
     56
     57   .. math::
     58
     59       K_\nu(x) = \frac{\pi}{2} \cdot \frac{I_{-\nu}(x) -
     60       I_n(x)}{\sin(\pi \nu)}.
     61
     62   We should note here that the above formulas for Bessel Y and K functions
     63   should be understood as limits when `\nu` is an integer.
     64
     65-  It follows from Bessel's differential equation that the derivative of
     66   `J_n(x)` with respect to `x` is:
     67
     68   .. math::
     69
     70       \frac{d}{dx} J_n(x) = \frac{1}{x^n} \left(x^n J_{n-1}(x) - n x^{n-1}
     71       J_n(z) \right)
     72
     73-  Another important formulation of the two linearly independent
     74   solutions to Bessel's equation are the Hankel functions
     75   `H_\nu^{(1)}(x)` and `H_\nu^{(2)}(x)`,
     76   defined by:
     77
     78   .. math::
     79
     80         H_\nu^{(1)}(x) = J_\nu(x) + i Y_\nu(x)
     81
     82   .. math::
     83
     84         H_\nu^{(2)}(x) = J_\nu(x) - i Y_\nu(x)
     85
     86   where `i` is the imaginary unit (and `J_*` and
     87   `Y_*` are the usual J- and Y-Bessel functions). These
     88   linear combinations are also known as Bessel functions of the third
     89   kind; they are also two linearly independent solutions of Bessel's
     90   differential equation. They are named for Hermann Hankel.
     91
     92EXAMPLES:
     93
     94    Evaluate the Bessel J function symbolically and numerically::
     95
     96        sage: bessel_J(0, x)
     97        bessel_J(0, x)
     98        sage: bessel_J(0, 0)
     99        bessel_J(0, 0)
     100        sage: bessel_J(0, x).diff(x)
     101        1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x)
     102
     103        sage: N(bessel_J(0, 0), digits = 20)
     104        1.0000000000000000000
     105        sage: find_root(bessel_J(0,x), 0, 5)
     106        2.404825557695773
     107
     108    Plot the Bessel J function::
     109
     110        sage: f(x) = Bessel(0)(x); f
     111        x |--> bessel_J(0, x)
     112        sage: plot(f, (x, 1, 10))
     113
     114    Visualize the Bessel Y function on the complex plane::
     115
     116        sage: complex_plot(bessel_Y(0, x), (-5, 5), (-5, 5))
     117
     118    Evaluate a combination of Bessel functions::
     119
     120        sage: f(x) = bessel_J(1, x) - bessel_Y(0, x)
     121        sage: f(pi)
     122        bessel_J(1, pi) - bessel_Y(0, pi)
     123        sage: f(pi).n()
     124        -0.0437509653365599
     125        sage: f(pi).n(digits=50)
     126        -0.043750965336559909054985168023342675387737118378169
     127
     128    Symbolically solve a second order differential equation with initial
     129    conditions `y(1) = a` and `y'(1) = b` in terms of Bessel functions::
     130
     131        sage: y = function('y', x)
     132        sage: a, b = var('a, b')
     133        sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0
     134        sage: f = desolve(diffeq, y, [1, a, b]); f
     135        (a*bessel_Y(1, 1) + b*bessel_Y(0, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (a*bessel_J(1, 1) + b*bessel_J(0, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))
     136
     137    For more examples, see the docstring for :meth:`Bessel`.
     138
     139AUTHORS:
     140
     141    - Benjamin Jones (2012-12-27): initial version.
     142
     143    - Some of the documentation here has been adapted from David Joyner's
     144      original documentation of Sage's special functions module (2006).
     145
     146REFERENCES:
     147
     148    - Abramowitz and Stegun: Handbook of Mathematical Functions,
     149      http://www.math.sfu.ca/~cbm/aands/
     150
     151    - http://en.wikipedia.org/wiki/Bessel_function
     152
     153    - mpmath Library `Bessel Functions`_
     154
     155.. _`mpmath Library`: http://code.google.com/p/mpmath/
     156.. _`Bessel Functions`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html
     157
     158"""
     159
     160#*****************************************************************************
     161#       Copyright (C) 2013 Benjamin Jones <benjaminfjones@gmail.com>
     162#
     163#  Distributed under the terms of the GNU General Public License (GPL)
     164#
     165#    This code is distributed in the hope that it will be useful,
     166#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     167#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     168#    General Public License for more details.
     169#
     170#  The full text of the GPL is available at:
     171#
     172#                  http://www.gnu.org/licenses/
     173#*****************************************************************************
     174
     175from sage.functions.other import sqrt
     176from sage.functions.log import exp
     177from sage.functions.hyperbolic import sinh, cosh
     178from sage.libs.mpmath import utils as mpmath_utils
     179from sage.misc.latex import latex
     180from sage.rings.all import RR, Integer
     181from sage.structure.coerce import parent
     182from sage.structure.element import get_coercion_model
     183from sage.symbolic.constants import pi
     184from sage.symbolic.function import BuiltinFunction, is_inexact
     185from sage.symbolic.expression import Expression
     186
     187# remove after deprecation period
     188from sage.calculus.calculus import maxima
     189from sage.functions.other import real, imag
     190from sage.misc.sage_eval import sage_eval
     191from sage.rings.real_mpfr import RealField
     192from sage.plot.plot import plot
     193from sage.rings.all import ZZ
     194
     195
     196class Function_Bessel_J(BuiltinFunction):
     197    r"""
     198    The Bessel J Function, denoted by bessel_J(`\nu`, x) or `J_\nu(x)`.
     199    As a Taylor series about `x=0` it is equal to:
     200
     201    .. math::
     202
     203        J_\nu(x) = \sum_{k=0}^\infty \frac{(-1)^k}{k! \Gamma(k+\nu+1)}
     204        (\frac{x}{2})^{2k+\nu}
     205
     206    The parameter `\nu` is called the order and may be any real or
     207    complex number, however the integer and half-integer values are most
     208    common. It is defined for all complex numbers `x` when `\nu`
     209    is an integer or greater than zero and it diverges as `x \to 0`
     210    for negative non-integer values of `\nu`.
     211
     212    For integer orders `\nu = n` there is an integral representation:
     213
     214    .. math::
     215
     216        J_n(x) = \frac{1}{\pi} \int_0^\pi \cos(n t - x \sin(t)) \; dt
     217
     218    This function also arises as a special case of the hypergeometric
     219    function `{}_0F_1`:
     220
     221    .. math::
     222
     223        J_\nu(x) = \frac{x^n}{2^\nu \Gamma(\nu + 1)} {}_0F_1(\nu +
     224        1, -\frac{x^2}{4}).
     225
     226    EXAMPLES::
     227
     228        sage: bessel_J(1.0, 1.0)
     229        0.440050585744933
     230        sage: bessel_J(2, I).n(digits=30)
     231        -0.135747669767038281182852569995
     232
     233        sage: bessel_J(1, x)
     234        bessel_J(1, x)
     235        sage: n = var('n')
     236        sage: bessel_J(n, x)
     237        bessel_J(n, x)
     238
     239    Examples of symbolic manipulation::
     240
     241        sage: a = bessel_J(pi, bessel_J(1, I)); a
     242        bessel_J(pi, bessel_J(1, I))
     243        sage: N(a, digits=20)
     244        0.00059023706363796717363 - 0.0026098820470081958110*I
     245
     246        sage: f = bessel_J(2, x)
     247        sage: f.diff(x)
     248        1/2*bessel_J(1, x) - 1/2*bessel_J(3, x)
     249
     250    Comparison to a well-known integral representation of `J_1(1)`::
     251
     252        sage: A = numerical_integral(1/pi*cos(x - sin(x)), 0, pi)
     253        sage: A[0]
     254        0.44005058574493355
     255        sage: bessel_J(1.0, 1.0) - A[0] < 1e-15
     256        True
     257
     258    Currently, integration is not supported (directly) since we cannot
     259    yet convert hypergeometric functions to and from Maxima::
     260
     261        sage: f = bessel_J(2, x)
     262        sage: f.integrate(x)
     263        Traceback (most recent call last):
     264        ...
     265        TypeError: cannot coerce arguments: no canonical coercion from <type 'list'> to Symbolic Ring
     266
     267        sage: m = maxima(bessel_J(2, x))
     268        sage: m.integrate(x)
     269        hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24
     270
     271    Visualization::
     272
     273        sage: plot(bessel_J(1,x), (x,0,5), color='blue')
     274        sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5)) # long time
     275
     276    ALGORITHM:
     277
     278        Numerical evaluation is handled by the mpmath library. Symbolics are
     279        handled by a combination of Maxima and Sage (Ginac/Pynac).
     280    """
     281    def __init__(self):
     282        """
     283        See the docstring for :meth:`Function_Bessel_J`.
     284
     285        EXAMPLES::
     286
     287            sage: sage.functions.bessel.Function_Bessel_J()
     288            bessel_J
     289        """
     290        BuiltinFunction.__init__(self, "bessel_J", nargs=2,
     291                                 conversions=dict(maxima='bessel_j',
     292                                                  mathematica='BesselJ'))
     293
     294    # remove after deprecation period
     295    def __call__(self, *args, **kwds):
     296        """
     297        Custom `__call__` method which uses the old Bessel function code if the
     298        `algorithm` or `prec` arguments are used. This should be removed after
     299        the deprecation period.
     300
     301        EXAMPLES::
     302
     303            sage: bessel_J(0, 1.0, "maxima", 53)
     304            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
     305            See http://trac.sagemath.org/4102 for details.
     306            .7651976865579666
     307        """
     308        if len(args) > 2 or len(kwds) > 0:
     309            from sage.misc.superseded import deprecation
     310            deprecation(4102, 'precision argument is deprecated; algorithm '
     311                              'argument is currently deprecated, but will be '
     312                              'available as a named keyword in the future')
     313            return _bessel_J(*args, **kwds)
     314        else:
     315            return super(BuiltinFunction, self).__call__(*args, **kwds)
     316
     317    def _eval_(self, n, x):
     318        """
     319        EXAMPLES::
     320
     321            sage: a, b = var('a, b')
     322            sage: bessel_J(a, b)
     323            bessel_J(a, b)
     324            sage: bessel_J(1.0, 1.0)
     325            0.440050585744933
     326        """
     327        if (not isinstance(n, Expression) and
     328                not isinstance(x, Expression) and
     329                (is_inexact(n) or is_inexact(x))):
     330            coercion_model = get_coercion_model()
     331            n, x = coercion_model.canonical_coercion(n, x)
     332            return self._evalf_(n, x, parent(n))
     333
     334        return None
     335
     336    def _evalf_(self, n, x, parent=None):
     337        """
     338        EXAMPLES::
     339
     340            sage: bessel_J(0.0, 1.0)
     341            0.765197686557966
     342            sage: bessel_J(0, 1).n(digits=20)
     343            0.76519768655796655145
     344        """
     345        import mpmath
     346        return mpmath_utils.call(mpmath.besselj, n, x, parent=parent)
     347
     348    def _derivative_(self, n, x, diff_param=None):
     349        """
     350        Return the derivative of the Bessel J function.
     351
     352        EXAMPLES::
     353
     354            sage: f(z) = bessel_J(10, z)
     355            sage: derivative(f, z)
     356            z |--> 1/2*bessel_J(9, z) - 1/2*bessel_J(11, z)
     357        """
     358        return (bessel_J(n-1, x) - bessel_J(n+1, x))/Integer(2)
     359
     360    def _print_latex_(self, n, z):
     361        """
     362        Custom _print_latex_ method.
     363
     364        EXAMPLES::
     365
     366            sage: latex(bessel_J(1, x))
     367            \operatorname{J_{1}}(x)
     368        """
     369        return r"\operatorname{J_{%s}}(%s)" % (latex(n), latex(z))
     370
     371bessel_J = Function_Bessel_J()
     372
     373
     374class Function_Bessel_Y(BuiltinFunction):
     375    r"""
     376    The Bessel Y functions, also known as the Bessel functions of the second
     377    kind, Weber functions, or Neumann functions.
     378
     379    `Y_\nu(z)` is a holomorphic function of `z` on the complex plane,
     380    cut along the negative real axis. It is singular at `z = 0`. When `z`
     381    is fixed, `Y_\nu(z)` is an entire function of the order `\nu`.
     382
     383    DEFINITION:
     384
     385    .. math::
     386
     387        Y_n(z) = \frac{J_\nu(z) \cos(\nu z) -
     388        J_{-\nu}(z)}{\sin(\nu z)}
     389
     390    Its derivative with respect to `z` is:
     391
     392    .. math::
     393
     394        \frac{d}{dz} Y_n(z) = \frac{1}{z^n} \left(z^n Y_{n-1}(z) - n z^{n-1}
     395        Y_n(z) \right)
     396
     397    EXAMPLES::
     398
     399        sage: bessel_Y(1, x)
     400        bessel_Y(1, x)
     401        sage: bessel_Y(1.0, 1.0)
     402        -0.781212821300289
     403        sage: n = var('n')
     404        sage: bessel_Y(n, x)
     405        bessel_Y(n, x)
     406        sage: bessel_Y(2, I).n()
     407        1.03440456978312 - 0.135747669767038*I
     408        sage: bessel_Y(0, 0).n()
     409        -infinity
     410
     411    Examples of symbolic manipulation::
     412
     413        sage: a = bessel_Y(pi, bessel_Y(1, I)); a
     414        bessel_Y(pi, bessel_Y(1, I))
     415        sage: N(a, digits=20)
     416        4.2059146571791095708 + 21.307914215321993526*I
     417
     418        sage: f = bessel_Y(2, x)
     419        sage: f.diff(x)
     420        1/2*bessel_Y(1, x) - 1/2*bessel_Y(3, x)
     421
     422    High precision and complex valued inputs (see :trac:`4230`)::
     423
     424        sage: bessel_Y(0, 1).n(128)
     425        0.088256964215676957982926766023515162828
     426        sage: bessel_Y(0, RealField(200)(1))
     427        0.088256964215676957982926766023515162827817523090675546711044
     428        sage: bessel_Y(0, ComplexField(200)(0.5+I))
     429        0.077763160184438051408593468823822434235010300228009867784073 + 1.0142336049916069152644677682828326441579314239591288411739*I
     430
     431    Visualization::
     432
     433        sage: plot(bessel_Y(1,x), (x,0,5), color='blue')
     434        sage: complex_plot(bessel_Y(1, x), (-5, 5), (-5, 5)) # long time
     435
     436    ALGORITHM:
     437
     438        Numerical evaluation is handled by the mpmath library. Symbolics are
     439        handled by a combination of Maxima and Sage (Ginac/Pynac).
     440    """
     441    def __init__(self):
     442        """
     443        See the docstring for :meth:`Function_Bessel_Y`.
     444
     445        EXAMPLES::
     446
     447            sage: sage.functions.bessel.Function_Bessel_Y()(0, x)
     448            bessel_Y(0, x)
     449        """
     450        BuiltinFunction.__init__(self, "bessel_Y", nargs=2,
     451                                 conversions=dict(maxima='bessel_y',
     452                                                  mathematica='BesselY'))
     453
     454    # remove after deprecation period
     455    def __call__(self, *args, **kwds):
     456        """
     457        Custom `__call__` method which uses the old Bessel function code if the
     458        `algorithm` or `prec` arguments are used. This should be removed after
     459        the deprecation period.
     460
     461        EXAMPLES::
     462
     463            sage: bessel_Y(0, 1, "maxima", 53)
     464            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
     465            See http://trac.sagemath.org/4102 for details.
     466            0.0882569642156769
     467        """
     468        if len(args) > 2 or len(kwds) > 0:
     469            from sage.misc.superseded import deprecation
     470            deprecation(4102, 'precision argument is deprecated; algorithm '
     471                              'argument is currently deprecated, but will be '
     472                              'available as a named keyword in the future')
     473            return _bessel_Y(*args, **kwds)
     474        else:
     475            return super(BuiltinFunction, self).__call__(*args, **kwds)
     476
     477    def _eval_(self, n, x):
     478        """
     479        EXAMPLES::
     480
     481            sage: a,b = var('a, b')
     482            sage: bessel_Y(a, b)
     483            bessel_Y(a, b)
     484            sage: bessel_Y(0, 1).n(128)
     485            0.088256964215676957982926766023515162828
     486        """
     487        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     488                (is_inexact(n) or is_inexact(x))):
     489            coercion_model = get_coercion_model()
     490            n, x = coercion_model.canonical_coercion(n, x)
     491            return self._evalf_(n, x, parent(n))
     492
     493        return None  # leaves the expression unevaluated
     494
     495    def _evalf_(self, n, x, parent=None):
     496        """
     497        EXAMPLES::
     498
     499            sage: bessel_Y(1.0+2*I, 3.0+4*I)
     500            0.699410324467538 + 0.228917940896421*I
     501            sage: bessel_Y(0, 1).n(256)
     502            0.08825696421567695798292676602351516282781752309067554671104384761199978932351
     503        """
     504        import mpmath
     505        return mpmath_utils.call(mpmath.bessely, n, x, parent=parent)
     506
     507    def _derivative_(self, n, x, diff_param=None):
     508        """
     509        Return the derivative of the Bessel Y function.
     510
     511        EXAMPLES::
     512
     513            sage: f(x) = bessel_Y(10, x)
     514            sage: derivative(f, x)
     515            x |--> 1/2*bessel_Y(9, x) - 1/2*bessel_Y(11, x)
     516        """
     517        return (bessel_Y(n-1, x) - bessel_Y(n+1, x))/Integer(2)
     518
     519    def _print_latex_(self, n, z):
     520        """
     521        Custom _print_latex_ method.
     522
     523        EXAMPLES::
     524
     525            sage: latex(bessel_Y(1, x))
     526            \operatorname{Y_{1}}(x)
     527        """
     528        return r"\operatorname{Y_{%s}}(%s)" % (latex(n), latex(z))
     529
     530bessel_Y = Function_Bessel_Y()
     531
     532
     533class Function_Bessel_I(BuiltinFunction):
     534    r"""
     535    The Bessel I function, or the Modified Bessel Function of the First Kind.
     536
     537    DEFINITION:
     538
     539    .. math::
     540
     541        I_\nu(x) = i^{-\nu} J_\nu(ix)
     542
     543    EXAMPLES::
     544
     545        sage: bessel_I(1, x)
     546        bessel_I(1, x)
     547        sage: bessel_I(1.0, 1.0)
     548        0.565159103992485
     549        sage: n = var('n')
     550        sage: bessel_I(n, x)
     551        bessel_I(n, x)
     552        sage: bessel_I(2, I).n()
     553        -0.114903484931900
     554
     555    Examples of symbolic manipulation::
     556
     557        sage: a = bessel_I(pi, bessel_I(1, I))
     558        sage: N(a, digits=20)
     559        0.00026073272117205890528 - 0.0011528954889080572266*I
     560
     561        sage: f = bessel_I(2, x)
     562        sage: f.diff(x)
     563        1/2*bessel_I(1, x) + 1/2*bessel_I(3, x)
     564
     565    Special identities that bessel_I satisfies::
     566
     567        sage: bessel_I(1/2, x)
     568        sqrt(1/(pi*x))*sqrt(2)*sinh(x)
     569        sage: eq = bessel_I(1/2, x) == bessel_I(0.5, x)
     570        sage: eq.test_relation()
     571        True
     572        sage: bessel_I(-1/2, x)
     573        sqrt(1/(pi*x))*sqrt(2)*cosh(x)
     574        sage: eq = bessel_I(-1/2, x) == bessel_I(-0.5, x)
     575        sage: eq.test_relation()
     576        True
     577
     578    Examples of asymptotic behavior::
     579
     580        sage: limit(bessel_I(0, x), x=oo)
     581        +Infinity
     582        sage: limit(bessel_I(0, x), x=0)
     583        1
     584
     585    High precision and complex valued inputs::
     586
     587        sage: bessel_I(0, 1).n(128)
     588        1.2660658777520083355982446252147175376
     589        sage: bessel_I(0, RealField(200)(1))
     590        1.2660658777520083355982446252147175376076703113549622068081
     591        sage: bessel_I(0, ComplexField(200)(0.5+I))
     592        0.80644357583493619472428518415019222845373366024179916785502 + 0.22686958987911161141397453401487525043310874687430711021434*I
     593
     594    Visualization::
     595
     596        sage: plot(bessel_I(1,x), (x,0,5), color='blue')
     597        sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5)) # long time
     598
     599    ALGORITHM:
     600
     601        Numerical evaluation is handled by the mpmath library. Symbolics are
     602        handled by a combination of Maxima and Sage (Ginac/Pynac).
     603    """
     604    def __init__(self):
     605        """
     606        See the docstring for :meth:`Function_Bessel_I`.
     607
     608        EXAMPLES::
     609
     610            sage: bessel_I(1,x)
     611            bessel_I(1, x)
     612        """
     613        BuiltinFunction.__init__(self, "bessel_I", nargs=2,
     614                                 conversions=dict(maxima='bessel_i',
     615                                                  mathematica='BesselI'))
     616
     617    # remove after deprecation period
     618    def __call__(self, *args, **kwds):
     619        """
     620        Custom `__call__` method which uses the old Bessel function code if the
     621        `algorithm` or `prec` arguments are used. This should be removed after
     622        the deprecation period.
     623
     624        EXAMPLES::
     625
     626            sage: bessel_I(0, 1, "maxima", 53)
     627            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
     628            See http://trac.sagemath.org/4102 for details.
     629            1.266065877752009
     630        """
     631        if len(args) > 2 or len(kwds) > 0:
     632            from sage.misc.superseded import deprecation
     633            deprecation(4102, 'precision argument is deprecated; algorithm '
     634                              'argument is currently deprecated, but will be '
     635                              'available as a named keyword in the future')
     636            return _bessel_I(*args, **kwds)
     637        else:
     638            return super(BuiltinFunction, self).__call__(*args, **kwds)
     639
     640    def _eval_(self, n, x):
     641        """
     642        EXAMPLES::
     643
     644            sage: y=var('y')
     645            sage: bessel_I(y,x)
     646            bessel_I(y, x)
     647            sage: bessel_I(0.0, 1.0)
     648            1.26606587775201
     649            sage: bessel_I(1/2, 1)
     650            sqrt(2)*sinh(1)/sqrt(pi)
     651            sage: bessel_I(-1/2, pi)
     652            sqrt(2)*cosh(pi)/pi
     653        """
     654        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     655                (is_inexact(n) or is_inexact(x))):
     656            coercion_model = get_coercion_model()
     657            n, x = coercion_model.canonical_coercion(n, x)
     658            return self._evalf_(n, x, parent(n))
     659
     660        # special identities
     661        if n == Integer(1)/Integer(2):
     662            return sqrt(2/(pi*x))*sinh(x)
     663        elif n == -Integer(1)/Integer(2):
     664            return sqrt(2/(pi*x))*cosh(x)
     665
     666        return None  # leaves the expression unevaluated
     667
     668    def _evalf_(self, n, x, parent=None):
     669        """
     670        EXAMPLES::
     671
     672            sage: bessel_I(1,3).n(digits=20)
     673            3.9533702174026093965
     674        """
     675        import mpmath
     676        return mpmath_utils.call(mpmath.besseli, n, x, parent=parent)
     677
     678    def _derivative_(self, n, x, diff_param=None):
     679        """
     680        Return the derivative of the Bessel I function `I_n(x)` with repect
     681        to `x`.
     682
     683        EXAMPLES::
     684
     685            sage: f(z) = bessel_I(10, x)
     686            sage: derivative(f, x)
     687            z |--> 1/2*bessel_I(9, x) + 1/2*bessel_I(11, x)
     688        """
     689        return (bessel_I(n-1, x) + bessel_I(n+1, x))/Integer(2)
     690
     691    def _print_latex_(self, n, z):
     692        """
     693        Custom _print_latex_ method.
     694
     695        EXAMPLES::
     696
     697            sage: latex(bessel_I(1, x))
     698            \operatorname{I_{1}}(x)
     699        """
     700        return r"\operatorname{I_{%s}}(%s)" % (latex(n), latex(z))
     701
     702bessel_I = Function_Bessel_I()
     703
     704
     705class Function_Bessel_K(BuiltinFunction):
     706    r"""
     707    The Bessel K function, or the modified Bessel function of the second kind.
     708
     709    DEFINITION:
     710
     711    .. math::
     712
     713        K_\nu(x) = \frac{\pi}{2} \frac{I_{-\nu}(x)-I_\nu(x)}{\sin(\nu \pi)}
     714
     715    EXAMPLES::
     716
     717        sage: bessel_K(1, x)
     718        bessel_K(1, x)
     719        sage: bessel_K(1.0, 1.0)
     720        0.601907230197235
     721        sage: n = var('n')
     722        sage: bessel_K(n, x)
     723        bessel_K(n, x)
     724        sage: bessel_K(2, I).n()
     725        -2.59288617549120 + 0.180489972066962*I
     726
     727    Examples of symbolic manipulation::
     728
     729        sage: a = bessel_K(pi, bessel_K(1, I)); a
     730        bessel_K(pi, bessel_K(1, I))
     731        sage: N(a, digits=20)
     732        3.8507583115005220157 + 0.068528298579883425792*I
     733
     734        sage: f = bessel_K(2, x)
     735        sage: f.diff(x)
     736        1/2*bessel_K(1, x) + 1/2*bessel_K(3, x)
     737
     738        sage: bessel_K(1/2, x)
     739        bessel_K(1/2, x)
     740        sage: bessel_K(1/2, -1)
     741        bessel_K(1/2, -1)
     742        sage: bessel_K(1/2, 1)
     743        sqrt(pi)*sqrt(1/2)*e^(-1)
     744
     745    Examples of asymptotic behavior::
     746
     747        sage: bessel_K(0, 0.0)
     748        +infinity
     749        sage: limit(bessel_K(0, x), x=0)
     750        +Infinity
     751        sage: limit(bessel_K(0, x), x=oo)
     752        0
     753
     754    High precision and complex valued inputs::
     755
     756        sage: bessel_K(0, 1).n(128)
     757        0.42102443824070833333562737921260903614
     758        sage: bessel_K(0, RealField(200)(1))
     759        0.42102443824070833333562737921260903613621974822666047229897
     760        sage: bessel_K(0, ComplexField(200)(0.5+I))
     761        0.058365979093103864080375311643360048144715516692187818271179 - 0.67645499731334483535184142196073004335768129348518210260256*I
     762
     763    Visualization::
     764
     765        sage: plot(bessel_K(1,x), (x,0,5), color='blue')
     766        sage: complex_plot(bessel_K(1, x), (-5, 5), (-5, 5)) # long time
     767
     768    ALGORITHM:
     769
     770        Numerical evaluation is handled by the mpmath library. Symbolics are
     771        handled by a combination of Maxima and Sage (Ginac/Pynac).
     772
     773    TESTS:
     774
     775    Verify that :trac:`3426` is fixed:
     776
     777    The Bessel K function can be evaluated numerically at complex orders::
     778
     779        sage: bessel_K(10 * I, 10).n()
     780        9.82415743819925e-8
     781
     782    For a fixed imaginary order and increasing, real, second component the
     783    value of Bessel K is exponentially decaying::
     784
     785        sage: for x in [10, 20, 50, 100, 200]: print bessel_K(5*I, x).n()
     786        5.27812176514912e-6
     787        3.11005908421801e-10
     788        2.66182488515423e-23 - 8.59622057747552e-58*I
     789        4.11189776828337e-45 - 1.01494840019482e-80*I
     790        1.15159692553603e-88 - 6.75787862113718e-125*I
     791    """
     792    def __init__(self):
     793        """
     794        See the docstring for :meth:`Function_Bessel_K`.
     795
     796        EXAMPLES::
     797
     798            sage: sage.functions.bessel.Function_Bessel_K()
     799            bessel_K
     800        """
     801        BuiltinFunction.__init__(self, "bessel_K", nargs=2,
     802                                 conversions=dict(maxima='bessel_k',
     803                                                  mathematica='BesselK'))
     804
     805    # remove after deprecation period
     806    def __call__(self, *args, **kwds):
     807        """
     808        Custom `__call__` method which uses the old Bessel function code if the
     809        `algorithm` or `prec` arguments are used. This should be removed after
     810        the deprecation period.
     811
     812        EXAMPLES::
     813
     814            sage: bessel_K(0, 1, "maxima", 53)
     815            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
     816            See http://trac.sagemath.org/4102 for details.
     817            0.0882569642156769
     818        """
     819        if len(args) > 2 or len(kwds) > 0:
     820            from sage.misc.superseded import deprecation
     821            deprecation(4102, 'precision argument is deprecated; algorithm '
     822                              'argument is currently deprecated, but will be '
     823                              'available as a named keyword in the future')
     824            return _bessel_Y(*args, **kwds)
     825        else:
     826            return super(BuiltinFunction, self).__call__(*args, **kwds)
     827
     828    def _eval_(self, n, x):
     829        """
     830        EXAMPLES::
     831
     832            sage: bessel_K(1,0)
     833            bessel_K(1, 0)
     834            sage: bessel_K(1.0, 0.0)
     835            +infinity
     836            sage: bessel_K(-1, 1).n(128)
     837            0.60190723019723457473754000153561733926
     838        """
     839        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     840                (is_inexact(n) or is_inexact(x))):
     841            coercion_model = get_coercion_model()
     842            n, x = coercion_model.canonical_coercion(n, x)
     843            return self._evalf_(n, x, parent(n))
     844
     845        # special identity
     846        if n == Integer(1)/Integer(2) and x > 0:
     847            return sqrt(pi/2)*exp(-x)*x**(Integer(1)/Integer(2))
     848
     849        return None  # leaves the expression unevaluated
     850
     851    def _evalf_(self, n, x, parent=None):
     852        """
     853        EXAMPLES::
     854
     855            sage: bessel_K(0.0, 1.0)
     856            0.421024438240708
     857            sage: bessel_K(0, RealField(128)(1))
     858            0.42102443824070833333562737921260903614
     859        """
     860        import mpmath
     861        return mpmath_utils.call(mpmath.besselk, n, x, parent=parent)
     862
     863    def _derivative_(self, n, x, diff_param=None):
     864        """
     865        Return the derivative of the Bessel K function.
     866
     867        EXAMPLES::
     868
     869            sage: f(x) = bessel_K(10, x)
     870            sage: derivative(f, x)
     871            x |--> 1/2*bessel_K(9, x) + 1/2*bessel_K(11, x)
     872        """
     873        return (bessel_K(n-1, x) + bessel_K(n+1, x))/Integer(2)
     874
     875    def _print_latex_(self, n, z):
     876        """
     877        Custom _print_latex_ method.
     878
     879        EXAMPLES::
     880
     881            sage: latex(bessel_K(1, x))
     882            \operatorname{K_{1}}(x)
     883        """
     884        return r"\operatorname{K_{%s}}(%s)" % (latex(n), latex(z))
     885
     886bessel_K = Function_Bessel_K()
     887
     888
     889# dictionary used in Bessel
     890bessel_type_dict = {'I': bessel_I, 'J': bessel_J, 'K': bessel_K, 'Y': bessel_Y}
     891
     892
     893def Bessel(*args, **kwds):
     894    """
     895    A function factory that produces symbolic I, J, K, and Y Bessel functions.
     896    There are several ways to call this function:
     897
     898        - ``Bessel(order, type)``
     899        - ``Bessel(order)`` -- type defaults to 'J'
     900        - ``Bessel(order, typ=T)``
     901        - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter
     902          function
     903        - ``Bessel()`` -- order is unspecified, type is 'J'
     904
     905    where `order` can be any integer and T must be one of the strings 'I', 'J',
     906    'K', or 'Y'.
     907
     908    See the EXAMPLES below.
     909
     910    EXAMPLES:
     911
     912    Construction of Bessel functions with various orders and types::
     913
     914        sage: Bessel()
     915        bessel_J
     916        sage: Bessel(1)(x)
     917        bessel_J(1, x)
     918        sage: Bessel(1, 'Y')(x)
     919        bessel_Y(1, x)
     920        sage: Bessel(-2, 'Y')(x)
     921        bessel_Y(-2, x)
     922        sage: Bessel(typ='K')
     923        bessel_K
     924        sage: Bessel(0, typ='I')(x)
     925        bessel_I(0, x)
     926
     927    Evaluation::
     928
     929        sage: f = Bessel(1)
     930        sage: f(3.0)
     931        0.339058958525936
     932        sage: f(3)
     933        bessel_J(1, 3)
     934        sage: f(3).n(digits=50)
     935        0.33905895852593645892551459720647889697308041819801
     936
     937        sage: g = Bessel(typ='J')
     938        sage: g(1,3)
     939        bessel_J(1, 3)
     940        sage: g(2, 3+I).n()
     941        0.634160370148554 + 0.0253384000032695*I
     942        sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15
     943        True
     944
     945    Symbolic calculus::
     946
     947        sage: f(x) = Bessel(0, 'J')(x)
     948        sage: derivative(f, x)
     949        x |--> 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x)
     950        sage: derivative(f, x, x)
     951        x |--> 1/4*bessel_J(-2, x) - 1/2*bessel_J(0, x) + 1/4*bessel_J(2, x)
     952
     953    Verify that `J_0` satisfies Bessel's differential equation numerically
     954    using the ``test_relation()`` method::
     955
     956        sage: y = bessel_J(0, x)
     957        sage: diffeq = x^2*derivative(y,x,x) + x*derivative(y,x) + x^2*y == 0
     958        sage: diffeq.test_relation(proof=False)
     959        True
     960
     961    Conversion to other systems::
     962
     963        sage: x,y = var('x,y')
     964        sage: f = maxima(Bessel(typ='K')(x,y))
     965        sage: f.derivative('x')
     966        %pi*csc(%pi*x)*('diff(bessel_i(-x,y),x,1)-'diff(bessel_i(x,y),x,1))/2-%pi*bessel_k(x,y)*cot(%pi*x)
     967        sage: f.derivative('y')
     968        -(bessel_k(x+1,y)+bessel_k(x-1,y))/2
     969
     970    Compute the particular solution to Bessel's Differential Equation that
     971    satisfies `y(1) = 1` and `y'(1) = 1`, then verify the initial conditions
     972    and plot it::
     973
     974        sage: y = function('y', x)
     975        sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0
     976        sage: f = desolve(diffeq, y, [1, 1, 1]); f
     977        (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))
     978
     979        sage: f.subs(x=1).n() # numerical verification
     980        1.00000000000000
     981        sage: fp = f.diff(x)
     982        sage: fp.subs(x=1).n()
     983        1.00000000000000
     984
     985        sage: f.subs(x=1).simplify_full() # symbolic verification
     986        1
     987        sage: fp = f.diff(x)
     988        sage: fp.subs(x=1).simplify_full()
     989        1
     990
     991        sage: plot(f, (x,0,5))
     992
     993    Plotting::
     994
     995        sage: f(x) = Bessel(0)(x); f
     996        x |--> bessel_J(0, x)
     997        sage: plot(f, (x, 1, 10))
     998
     999        sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10)
     1000
     1001        sage: G = Graphics()
     1002        sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ])
     1003        sage: show(G)
     1004
     1005    A recreation of Abramowitz and Stegun Figure 9.1::
     1006
     1007        sage: G  = plot(Bessel(0, 'J'), 0, 15, color='black')
     1008        sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black')
     1009        sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted')
     1010        sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted')
     1011        sage: show(G, ymin=-1, ymax=1)
     1012
     1013    """
     1014    # Determine the order and type of function from the arguments and keywords.
     1015    # These are recored in local variables: _type, _order, _system, _nargs.
     1016    _type = None
     1017    if len(args) == 0:    # no order specified
     1018        _order = None
     1019        _nargs = 2
     1020    elif len(args) == 1:  # order is specified
     1021        _order = args[0]
     1022        _nargs = 1
     1023    elif len(args) == 2:  # both order and type are positional arguments
     1024        _order = args[0]
     1025        _type = args[1]
     1026        _nargs = 1
     1027    else:
     1028        from sage.misc.superseded import deprecation
     1029        deprecation(4102, 'precision argument is deprecated; algorithm '
     1030                          'argument is currently deprecated, but will be '
     1031                          'available as a named keyword in the future')
     1032        return _Bessel(*args, **kwds)
     1033
     1034    # check for type inconsistency
     1035    if _type is not None and 'typ' in kwds and _type != kwds['typ']:
     1036        raise ValueError("inconsistent types given")
     1037    # record the function type
     1038    if _type is None:
     1039        if 'typ' in kwds:
     1040            _type = kwds['typ']
     1041        else:
     1042            _type = 'J'
     1043    if not (_type in ['I', 'J', 'K', 'Y']):
     1044        raise ValueError("type must be one of I, J, K, Y")
     1045    # record the numerical evaluation system
     1046    if 'algorithm' in kwds:
     1047        _system = kwds['algorithm']
     1048    else:
     1049        _system = 'mpmath'
     1050
     1051    # return the function
     1052    # TODO remove assertions
     1053    assert _type in ['I', 'J', 'K', 'Y']
     1054    assert _order is None or _order in RR
     1055    assert _nargs == 1 or _nargs == 2
     1056    assert _system == 'mpmath'
     1057
     1058    # TODO what to do with _system?
     1059    _f = bessel_type_dict[_type]
     1060    if _nargs == 1:
     1061        return lambda x: _f(_order, x)
     1062    else:
     1063        return _f
     1064
     1065####################################################
     1066###  to be removed after the deprecation period  ###
     1067####################################################
     1068
     1069
     1070def _bessel_I(nu,z,algorithm = "pari",prec=53):
     1071    r"""
     1072    Implements the "I-Bessel function", or "modified Bessel function,
     1073    1st kind", with index (or "order") nu and argument z.
     1074
     1075    INPUT:
     1076
     1077
     1078    -  ``nu`` - a real (or complex, for pari) number
     1079
     1080    -  ``z`` - a real (positive) algorithm - "pari" or
     1081       "maxima" or "scipy" prec - real precision (for PARI only)
     1082
     1083
     1084    DEFINITION::
     1085
     1086            Maxima:
     1087                             inf
     1088                            ====   - nu - 2 k  nu + 2 k
     1089                            \     2          z
     1090                             >    -------------------
     1091                            /     k! Gamma(nu + k + 1)
     1092                            ====
     1093                            k = 0
     1094
     1095            PARI:
     1096
     1097                             inf
     1098                            ====   - 2 k  2 k
     1099                            \     2      z    Gamma(nu + 1)
     1100                             >    -----------------------
     1101                            /       k! Gamma(nu + k + 1)
     1102                            ====
     1103                            k = 0
     1104
     1105
     1106
     1107    Sometimes ``bessel_I(nu,z)`` is denoted
     1108    ``I_nu(z)`` in the literature.
     1109
     1110    .. warning::
     1111
     1112       In Maxima (the manual says) i0 is deprecated but
     1113       ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
     1114       though.)
     1115
     1116    EXAMPLES::
     1117
     1118        sage: from sage.functions.bessel import _bessel_I
     1119        sage: _bessel_I(1,1,"pari",500)
     1120        0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
     1121        sage: _bessel_I(1,1)
     1122        0.565159103992485
     1123        sage: _bessel_I(2,1.1,"maxima")
     1124        0.16708949925104...
     1125        sage: _bessel_I(0,1.1,"maxima")
     1126        1.32616018371265...
     1127        sage: _bessel_I(0,1,"maxima")
     1128        1.2660658777520...
     1129        sage: _bessel_I(1,1,"scipy")
     1130        0.565159103992...
     1131
     1132    Check whether the return value is real whenever the argument is real (#10251)::
     1133
     1134        sage: _bessel_I(5, 1.5, algorithm='scipy') in RR
     1135        True
     1136
     1137    """
     1138    if algorithm=="pari":
     1139        from sage.libs.pari.all import pari
     1140        try:
     1141            R = RealField(prec)
     1142            nu = R(nu)
     1143            z = R(z)
     1144        except TypeError:
     1145            C = ComplexField(prec)
     1146            nu = C(nu)
     1147            z = C(z)
     1148            K = C
     1149        K = z.parent()
     1150        return K(pari(nu).besseli(z, precision=prec))
     1151    elif algorithm=="scipy":
     1152        if prec != 53:
     1153            raise ValueError, "for the scipy algorithm the precision must be 53"
     1154        import scipy.special
     1155        ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
     1156        ans = ans.replace("(","")
     1157        ans = ans.replace(")","")
     1158        ans = ans.replace("j","*I")
     1159        ans = sage_eval(ans)
     1160        return real(ans) if z in RR else ans # Return real value when arg is real
     1161    elif algorithm == "maxima":
     1162        if prec != 53:
     1163            raise ValueError, "for the maxima algorithm the precision must be 53"
     1164        return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
     1165    else:
     1166        raise ValueError, "unknown algorithm '%s'"%algorithm
     1167
     1168def _bessel_J(nu,z,algorithm="pari",prec=53):
     1169    r"""
     1170    Return value of the "J-Bessel function", or "Bessel function, 1st
     1171    kind", with index (or "order") nu and argument z.
     1172
     1173    ::
     1174
     1175            Defn:
     1176            Maxima:
     1177                             inf
     1178                            ====          - nu - 2 k  nu + 2 k
     1179                            \     (-1)^k 2           z
     1180                             >    -------------------------
     1181                            /        k! Gamma(nu + k + 1)
     1182                            ====
     1183                            k = 0
     1184
     1185            PARI:
     1186
     1187                             inf
     1188                            ====          - 2k    2k
     1189                            \     (-1)^k 2      z    Gamma(nu + 1)
     1190                             >    ----------------------------
     1191                            /         k! Gamma(nu + k + 1)
     1192                            ====
     1193                            k = 0
     1194
     1195
     1196    Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
     1197
     1198    .. warning::
     1199
     1200       Inaccurate for small values of z.
     1201
     1202    EXAMPLES::
     1203
     1204        sage: from sage.functions.bessel import _bessel_J
     1205        sage: _bessel_J(2,1.1)
     1206        0.136564153956658
     1207        sage: _bessel_J(0,1.1)
     1208        0.719622018527511
     1209        sage: _bessel_J(0,1)
     1210        0.765197686557967
     1211        sage: _bessel_J(0,0)
     1212        1.00000000000000
     1213        sage: _bessel_J(0.1,0.1)
     1214        0.777264368097005
     1215
     1216    We check consistency of PARI and Maxima::
     1217
     1218        sage: n(_bessel_J(3,10,"maxima"))
     1219        0.0583793793051...
     1220        sage: n(_bessel_J(3,10,"pari"))
     1221        0.0583793793051868
     1222        sage: _bessel_J(3,10,"scipy")
     1223        0.0583793793052...
     1224
     1225    Check whether the return value is real whenever the argument is real (#10251)::
     1226        sage: _bessel_J(5, 1.5, algorithm='scipy') in RR
     1227        True
     1228    """
     1229
     1230    if algorithm=="pari":
     1231        from sage.libs.pari.all import pari
     1232        try:
     1233            R = RealField(prec)
     1234            nu = R(nu)
     1235            z = R(z)
     1236        except TypeError:
     1237            C = ComplexField(prec)
     1238            nu = C(nu)
     1239            z = C(z)
     1240            K = C
     1241        if nu == 0:
     1242            nu = ZZ(0)
     1243        K = z.parent()
     1244        return K(pari(nu).besselj(z, precision=prec))
     1245    elif algorithm=="scipy":
     1246        if prec != 53:
     1247            raise ValueError, "for the scipy algorithm the precision must be 53"
     1248        import scipy.special
     1249        ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
     1250        ans = ans.replace("(","")
     1251        ans = ans.replace(")","")
     1252        ans = ans.replace("j","*I")
     1253        ans = sage_eval(ans)
     1254        return real(ans) if z in RR else ans
     1255    elif algorithm == "maxima":
     1256        if prec != 53:
     1257            raise ValueError, "for the maxima algorithm the precision must be 53"
     1258        f = maxima.function('n,z', 'bessel_j(n, z)')
     1259        return f(nu, z)
     1260    else:
     1261        raise ValueError, "unknown algorithm '%s'"%algorithm
     1262
     1263def _bessel_K(nu,z,algorithm="pari",prec=53):
     1264    r"""
     1265    Implements the "K-Bessel function", or "modified Bessel function,
     1266    2nd kind", with index (or "order") nu and argument z. Defn::
     1267
     1268                    pi*(bessel_I(-nu, z) - bessel_I(nu, z))
     1269                   ----------------------------------------
     1270                                2*sin(pi*nu)
     1271
     1272
     1273    if nu is not an integer and by taking a limit otherwise.
     1274
     1275    Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
     1276    PARI, nu can be complex and z must be real and positive.
     1277
     1278    EXAMPLES::
     1279
     1280        sage: from sage.functions.bessel import _bessel_K
     1281        sage: _bessel_K(3,2,"scipy")
     1282        0.64738539094...
     1283        sage: _bessel_K(3,2)
     1284        0.64738539094...
     1285        sage: _bessel_K(1,1)
     1286        0.60190723019...
     1287        sage: _bessel_K(1,1,"pari",10)
     1288        0.60
     1289        sage: _bessel_K(1,1,"pari",100)
     1290        0.60190723019723457473754000154
     1291
     1292    TESTS::
     1293
     1294        sage: _bessel_K(2,1.1, algorithm="maxima")
     1295        Traceback (most recent call last):
     1296        ...
     1297        NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms
     1298
     1299        Check whether the return value is real whenever the argument is real (#10251)::
     1300
     1301        sage: _bessel_K(5, 1.5, algorithm='scipy') in RR
     1302        True
     1303
     1304    """
     1305    if algorithm=="scipy":
     1306        if prec != 53:
     1307            raise ValueError, "for the scipy algorithm the precision must be 53"
     1308        import scipy.special
     1309        ans = str(scipy.special.kv(float(nu),float(z)))
     1310        ans = ans.replace("(","")
     1311        ans = ans.replace(")","")
     1312        ans = ans.replace("j","*I")
     1313        ans = sage_eval(ans)
     1314        return real(ans) if z in RR else ans
     1315    elif algorithm == 'pari':
     1316        from sage.libs.pari.all import pari
     1317        try:
     1318            R = RealField(prec)
     1319            nu = R(nu)
     1320            z = R(z)
     1321        except TypeError:
     1322            C = ComplexField(prec)
     1323            nu = C(nu)
     1324            z = C(z)
     1325            K = C
     1326        K = z.parent()
     1327        return K(pari(nu).besselk(z, precision=prec))
     1328    elif algorithm == 'maxima':
     1329        raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
     1330    else:
     1331        raise ValueError, "unknown algorithm '%s'"%algorithm
     1332
     1333
     1334def _bessel_Y(nu,z,algorithm="maxima", prec=53):
     1335    r"""
     1336    Implements the "Y-Bessel function", or "Bessel function of the 2nd
     1337    kind", with index (or "order") nu and argument z.
     1338
     1339    .. note::
     1340
     1341       Currently only prec=53 is supported.
     1342
     1343    Defn::
     1344
     1345                    cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
     1346                   -------------------------------------------------
     1347                                     sin(nu*pi)
     1348
     1349    if nu is not an integer and by taking a limit otherwise.
     1350
     1351    Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
     1352
     1353    This is computed using Maxima by default.
     1354
     1355    EXAMPLES::
     1356
     1357        sage: from sage.functions.bessel import _bessel_Y
     1358        sage: _bessel_Y(2,1.1,"scipy")
     1359        -1.4314714939...
     1360        sage: _bessel_Y(2,1.1)
     1361        -1.4314714939590...
     1362        sage: _bessel_Y(3.001,2.1)
     1363        -1.0299574976424...
     1364
     1365    TESTS::
     1366
     1367        sage: _bessel_Y(2,1.1, algorithm="pari")
     1368        Traceback (most recent call last):
     1369        ...
     1370        NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
     1371    """
     1372    if algorithm=="scipy":
     1373        if prec != 53:
     1374            raise ValueError, "for the scipy algorithm the precision must be 53"
     1375        import scipy.special
     1376        ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
     1377        ans = ans.replace("(","")
     1378        ans = ans.replace(")","")
     1379        ans = ans.replace("j","*I")
     1380        ans = sage_eval(ans)
     1381        return real(ans) if z in RR else ans
     1382    elif algorithm == "maxima":
     1383        if prec != 53:
     1384            raise ValueError, "for the maxima algorithm the precision must be 53"
     1385        return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
     1386    elif algorithm == "pari":
     1387        raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
     1388    else:
     1389        raise ValueError, "unknown algorithm '%s'"%algorithm
     1390
     1391class _Bessel():
     1392    """
     1393    A class implementing the I, J, K, and Y Bessel functions.
     1394
     1395    EXAMPLES::
     1396
     1397        sage: from sage.functions.bessel import _Bessel
     1398        sage: g = _Bessel(2); g
     1399        J_{2}
     1400        sage: print g
     1401        J-Bessel function of order 2
     1402        sage: g.plot(0,10)
     1403
     1404    ::
     1405
     1406        sage: _Bessel(2, typ='I')(pi)
     1407        2.61849485263445
     1408        sage: _Bessel(2, typ='J')(pi)
     1409        0.485433932631509
     1410        sage: _Bessel(2, typ='K')(pi)
     1411        0.0510986902537926
     1412        sage: _Bessel(2, typ='Y')(pi)
     1413        -0.0999007139289404
     1414    """
     1415    def __init__(self, nu, typ = "J", algorithm = None, prec = 53):
     1416        """
     1417        Initializes new instance of the Bessel class.
     1418
     1419        INPUT:
     1420
     1421         - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'
     1422           or 'Y'.
     1423
     1424         - ``algorithm`` -- (default: maxima for type Y, pari for other types)
     1425           algorithm to use to compute the Bessel function: 'pari', 'maxima' or
     1426           'scipy'.  Note that type K is not implemented in Maxima and type Y
     1427           is not implemented in PARI.
     1428
     1429         - ``prec`` -- (default: 53) precision in bits of the Bessel function.
     1430           Only supported for the PARI algorithm.
     1431
     1432        EXAMPLES::
     1433
     1434            sage: from sage.functions.bessel import _Bessel
     1435            sage: g = _Bessel(2); g
     1436            J_{2}
     1437            sage: _Bessel(1,'I')
     1438            I_{1}
     1439            sage: _Bessel(6, prec=120)(pi)
     1440            0.014545966982505560573660369604001804
     1441            sage: _Bessel(6, algorithm="pari")(pi)
     1442            0.0145459669825056
     1443
     1444        For the Bessel J-function, Maxima returns a symbolic result.  For
     1445        types I and Y, we always get a numeric result::
     1446
     1447            sage: b = _Bessel(6, algorithm="maxima")(pi); b
     1448            bessel_j(6,pi)
     1449            sage: b.n(53)
     1450            0.0145459669825056
     1451            sage: _Bessel(6, typ='I', algorithm="maxima")(pi)
     1452            0.0294619840059568
     1453            sage: _Bessel(6, typ='Y', algorithm="maxima")(pi)
     1454            -4.33932818939038
     1455
     1456        SciPy usually gives less precise results::
     1457
     1458            sage: _Bessel(6, algorithm="scipy")(pi)
     1459            0.0145459669825000...
     1460
     1461        TESTS::
     1462
     1463            sage: _Bessel(1,'Z')
     1464            Traceback (most recent call last):
     1465            ...
     1466            ValueError: typ must be one of I, J, K, Y
     1467        """
     1468        if not (typ in ['I', 'J', 'K', 'Y']):
     1469            raise ValueError, "typ must be one of I, J, K, Y"
     1470
     1471        # Did the user ask for the default algorithm?
     1472        if algorithm is None:
     1473            if typ == 'Y':
     1474                algorithm = 'maxima'
     1475            else:
     1476                algorithm = 'pari'
     1477
     1478        self._system = algorithm
     1479        self._order = nu
     1480        self._type = typ
     1481        prec = int(prec)
     1482        if prec < 0:
     1483            raise ValueError, "prec must be a positive integer"
     1484        self._prec = int(prec)
     1485
     1486    def __str__(self):
     1487        """
     1488        Returns a string representation of this Bessel object.
     1489
     1490        TEST::
     1491
     1492            sage: from sage.functions.bessel import _Bessel
     1493            sage: a = _Bessel(1,'I')
     1494            sage: str(a)
     1495            'I-Bessel function of order 1'
     1496        """
     1497        return self.type()+"-Bessel function of order "+str(self.order())
     1498
     1499    def __repr__(self):
     1500        """
     1501        Returns a string representation of this Bessel object.
     1502
     1503        TESTS::
     1504
     1505            sage: from sage.functions.bessel import _Bessel
     1506            sage: _Bessel(1,'I')
     1507            I_{1}
     1508        """
     1509        return self.type()+"_{"+str(self.order())+"}"
     1510
     1511    def type(self):
     1512        """
     1513        Returns the type of this Bessel object.
     1514
     1515        TEST::
     1516
     1517            sage: from sage.functions.bessel import _Bessel
     1518            sage: a = _Bessel(3,'K')
     1519            sage: a.type()
     1520            'K'
     1521        """
     1522        return self._type
     1523
     1524    def prec(self):
     1525        """
     1526        Returns the precision (in number of bits) used to represent this
     1527        Bessel function.
     1528
     1529        TESTS::
     1530
     1531            sage: from sage.functions.bessel import _Bessel
     1532            sage: a = _Bessel(3,'K')
     1533            sage: a.prec()
     1534            53
     1535            sage: B = _Bessel(20,prec=100); B
     1536            J_{20}
     1537            sage: B.prec()
     1538            100
     1539        """
     1540        return self._prec
     1541
     1542    def order(self):
     1543        """
     1544        Returns the order of this Bessel function.
     1545
     1546        TEST::
     1547
     1548            sage: from sage.functions.bessel import _Bessel
     1549            sage: a = _Bessel(3,'K')
     1550            sage: a.order()
     1551            3
     1552        """
     1553        return self._order
     1554
     1555    def system(self):
     1556        """
     1557        Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with
     1558        this Bessel function.
     1559
     1560        TESTS::
     1561
     1562            sage: from sage.functions.bessel import _Bessel
     1563            sage: _Bessel(20,algorithm='maxima').system()
     1564            'maxima'
     1565            sage: _Bessel(20,prec=100).system()
     1566            'pari'
     1567        """
     1568        return self._system
     1569
     1570    def __call__(self,z):
     1571        """
     1572        Implements evaluation of all the Bessel functions directly
     1573        from the Bessel class. This essentially allows a Bessel object to
     1574        behave like a function that can be invoked.
     1575
     1576        TESTS::
     1577
     1578            sage: from sage.functions.bessel import _Bessel
     1579            sage: _Bessel(3,'K')(5.0)
     1580            0.00829176841523093
     1581            sage: _Bessel(20,algorithm='maxima')(5.0)
     1582            27.703300521289436e-12
     1583            sage: _Bessel(20,prec=100)(5.0101010101010101)
     1584            2.8809188227195382093062257967e-11
     1585            sage: B = _Bessel(2,'Y',algorithm='scipy',prec=50)
     1586            sage: B(2.0)
     1587            Traceback (most recent call last):
     1588            ...
     1589            ValueError: for the scipy algorithm the precision must be 53
     1590        """
     1591        nu = self.order()
     1592        t = self.type()
     1593        s = self.system()
     1594        p = self.prec()
     1595        if t == "I":
     1596            return _bessel_I(nu,z,algorithm=s,prec=p)
     1597        if t == "J":
     1598            return _bessel_J(nu,z,algorithm=s,prec=p)
     1599        if t == "K":
     1600            return _bessel_K(nu,z,algorithm=s,prec=p)
     1601        if t == "Y":
     1602            return _bessel_Y(nu,z,algorithm=s,prec=p)
     1603
     1604    def plot(self,a,b):
     1605        """
     1606        Enables easy plotting of all the Bessel functions directly
     1607        from the Bessel class.
     1608
     1609        TESTS::
     1610
     1611            sage: from sage.functions.bessel import _Bessel
     1612            sage: plot(_Bessel(2),3,4)
     1613            sage: _Bessel(2).plot(3,4)
     1614            sage: P = _Bessel(2,'I').plot(1,5)
     1615            sage: P += _Bessel(2,'J').plot(1,5)
     1616            sage: P += _Bessel(2,'K').plot(1,5)
     1617            sage: P += _Bessel(2,'Y').plot(1,5)
     1618            sage: show(P)
     1619        """
     1620        nu = self.order()
     1621        s = self.system()
     1622        t = self.type()
     1623        if t == "I":
     1624            f = lambda z: _bessel_I(nu,z,s)
     1625            P = plot(f,a,b)
     1626        if t == "J":
     1627            f = lambda z: _bessel_J(nu,z,s)
     1628            P = plot(f,a,b)
     1629        if t == "K":
     1630            f = lambda z: _bessel_K(nu,z,s)
     1631            P = plot(f,a,b)
     1632        if t == "Y":
     1633            f = lambda z: _bessel_Y(nu,z,s)
     1634            P = plot(f,a,b)
     1635        return P
  • sage/functions/special.py

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    2626Toy. It is placed under the terms of the General Public License
    2727(GPL) that governs the distribution of Maxima.
    2828
    29 The (usual) Bessel functions and Airy functions are part of the
    30 standard Maxima package. Some Bessel functions also are implemented
    31 in PARI. (Caution: The PARI versions are sometimes different than
    32 the Maxima version.) For example, the J-Bessel function
    33 `J_\nu (z)` can be computed using either Maxima or PARI,
    34 depending on an optional variable you pass to bessel_J.
    35 
    3629Next, we summarize some of the properties of the functions
    3730implemented here.
    3831
    3932
    40 -  Bessel functions, first defined by the Swiss mathematician
    41    Daniel Bernoulli and named after Friedrich Bessel, are canonical
    42    solutions y(x) of Bessel's differential equation:
    43 
    44    
    45    .. math::
    46 
    47          x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0,
    48 
    49 
    50    for an arbitrary real number `\alpha` (the order).
    51 
    52 -  Another important formulation of the two linearly independent
    53    solutions to Bessel's equation are the Hankel functions
    54    `H_\alpha^{(1)}(x)` and `H_\alpha^{(2)}(x)`,
    55    defined by:
    56 
    57    
    58    .. math::
    59 
    60          H_\alpha^{(1)}(x) = J_\alpha(x) + i Y_\alpha(x)
    61 
    62 
    63    
    64    .. math::
    65 
    66          H_\alpha^{(2)}(x) = J_\alpha(x) - i Y_\alpha(x)
    67 
    68 
    69    where `i` is the imaginary unit (and `J_*` and
    70    `Y_*` are the usual J- and Y-Bessel functions). These
    71    linear combinations are also known as Bessel functions of the third
    72    kind; they are two linearly independent solutions of Bessel's
    73    differential equation. They are named for Hermann Hankel.
    74 
    7533-  Airy function The function `Ai(x)` and the related
    7634   function `Bi(x)`, which is also called an Airy function,
    7735   are solutions to the differential equation
     
    333291- Abramowitz and Stegun: Handbook of Mathematical Functions,
    334292  http://www.math.sfu.ca/~cbm/aands/
    335293
    336 - http://en.wikipedia.org/wiki/Bessel_function
    337 
    338294- http://en.wikipedia.org/wiki/Airy_function
    339295
    340296- http://en.wikipedia.org/wiki/Spherical_harmonics
     
    567523
    568524    EXAMPLES::
    569525
    570         sage: n(bessel_J(3,10,"maxima"))
    571         0.0583793793051...
    572526        sage: spherical_hankel2(2,i)
    573527        -e
    574528    """
     
    582536
    583537            TESTS::
    584538
    585                 sage: n(bessel_J(3,10,"maxima"))
    586                 0.0583793793051...
    587539                sage: spherical_hankel2(2,x)
    588540                (-I*x^2 - 3*x + 3*I)*e^(-I*x)/x^3
    589541            """
     
    665617   return RDF(meval("airy_bi(%s)"%RDF(x)))
    666618
    667619
    668 def bessel_I(nu,z,algorithm = "pari",prec=53):
    669     r"""
    670     Implements the "I-Bessel function", or "modified Bessel function,
    671     1st kind", with index (or "order") nu and argument z.
    672    
    673     INPUT:
    674    
    675    
    676     -  ``nu`` - a real (or complex, for pari) number
    677    
    678     -  ``z`` - a real (positive) algorithm - "pari" or
    679        "maxima" or "scipy" prec - real precision (for PARI only)
    680    
    681    
    682     DEFINITION::
    683    
    684             Maxima:
    685                              inf
    686                             ====   - nu - 2 k  nu + 2 k
    687                             \     2          z
    688                              >    -------------------
    689                             /     k! Gamma(nu + k + 1)
    690                             ====
    691                             k = 0
    692        
    693             PARI:
    694            
    695                              inf
    696                             ====   - 2 k  2 k
    697                             \     2      z    Gamma(nu + 1)
    698                              >    -----------------------
    699                             /       k! Gamma(nu + k + 1)
    700                             ====
    701                             k = 0
    702        
    703            
    704    
    705     Sometimes ``bessel_I(nu,z)`` is denoted
    706     ``I_nu(z)`` in the literature.
    707    
    708     .. warning::
    709 
    710        In Maxima (the manual says) i0 is deprecated but
    711        ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
    712        though.)
    713    
    714     EXAMPLES::
    715    
    716         sage: bessel_I(1,1,"pari",500)
    717         0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
    718         sage: bessel_I(1,1)
    719         0.565159103992485
    720         sage: bessel_I(2,1.1,"maxima") 
    721         0.16708949925104...
    722         sage: bessel_I(0,1.1,"maxima")
    723         1.32616018371265...
    724         sage: bessel_I(0,1,"maxima")   
    725         1.2660658777520...
    726         sage: bessel_I(1,1,"scipy")
    727         0.565159103992...
    728 
    729     Check whether the return value is real whenever the argument is real (#10251)::
    730    
    731         sage: bessel_I(5, 1.5, algorithm='scipy') in RR
    732         True
    733        
    734     """
    735     if algorithm=="pari":
    736         from sage.libs.pari.all import pari
    737         try:
    738             R = RealField(prec)
    739             nu = R(nu)
    740             z = R(z)
    741         except TypeError:
    742             C = ComplexField(prec)
    743             nu = C(nu)
    744             z = C(z)
    745             K = C
    746         K = z.parent()
    747         return K(pari(nu).besseli(z, precision=prec))
    748     elif algorithm=="scipy":
    749         if prec != 53:
    750             raise ValueError, "for the scipy algorithm the precision must be 53"
    751         import scipy.special
    752         ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
    753         ans = ans.replace("(","")
    754         ans = ans.replace(")","")
    755         ans = ans.replace("j","*I")
    756         ans = sage_eval(ans)
    757         return real(ans) if z in RR else ans # Return real value when arg is real
    758     elif algorithm == "maxima":
    759         if prec != 53:
    760             raise ValueError, "for the maxima algorithm the precision must be 53"
    761         return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
    762     else:
    763         raise ValueError, "unknown algorithm '%s'"%algorithm
    764        
    765 def bessel_J(nu,z,algorithm="pari",prec=53):
    766     r"""
    767     Return value of the "J-Bessel function", or "Bessel function, 1st
    768     kind", with index (or "order") nu and argument z.
    769    
    770     ::
    771    
    772             Defn:
    773             Maxima:
    774                              inf
    775                             ====          - nu - 2 k  nu + 2 k
    776                             \     (-1)^k 2           z
    777                              >    -------------------------
    778                             /        k! Gamma(nu + k + 1)
    779                             ====
    780                             k = 0
    781        
    782             PARI:
    783            
    784                              inf
    785                             ====          - 2k    2k
    786                             \     (-1)^k 2      z    Gamma(nu + 1)
    787                              >    ----------------------------
    788                             /         k! Gamma(nu + k + 1)
    789                             ====
    790                             k = 0
    791            
    792    
    793     Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
    794    
    795     .. warning::
    796 
    797        Inaccurate for small values of z.
    798    
    799     EXAMPLES::
    800    
    801         sage: bessel_J(2,1.1)
    802         0.136564153956658
    803         sage: bessel_J(0,1.1)
    804         0.719622018527511
    805         sage: bessel_J(0,1)
    806         0.765197686557967
    807         sage: bessel_J(0,0)
    808         1.00000000000000
    809         sage: bessel_J(0.1,0.1)
    810         0.777264368097005
    811    
    812     We check consistency of PARI and Maxima::
    813    
    814         sage: n(bessel_J(3,10,"maxima"))
    815         0.0583793793051...
    816         sage: n(bessel_J(3,10,"pari")) 
    817         0.0583793793051868
    818         sage: bessel_J(3,10,"scipy")
    819         0.0583793793052...
    820 
    821     Check whether the return value is real whenever the argument is real (#10251)::                                                                                                                                                           
    822         sage: bessel_J(5, 1.5, algorithm='scipy') in RR                                                                     
    823         True
    824     """
    825    
    826     if algorithm=="pari":
    827         from sage.libs.pari.all import pari
    828         try:
    829             R = RealField(prec)
    830             nu = R(nu)
    831             z = R(z)
    832         except TypeError:
    833             C = ComplexField(prec)
    834             nu = C(nu)
    835             z = C(z)
    836             K = C
    837         if nu == 0:
    838             nu = ZZ(0)
    839         K = z.parent()
    840         return K(pari(nu).besselj(z, precision=prec))
    841     elif algorithm=="scipy":
    842         if prec != 53:
    843             raise ValueError, "for the scipy algorithm the precision must be 53"
    844         import scipy.special
    845         ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
    846         ans = ans.replace("(","")
    847         ans = ans.replace(")","")
    848         ans = ans.replace("j","*I")
    849         ans = sage_eval(ans)
    850         return real(ans) if z in RR else ans
    851     elif algorithm == "maxima":
    852         if prec != 53:
    853             raise ValueError, "for the maxima algorithm the precision must be 53"
    854         return maxima_function("bessel_j")(nu, z)
    855     else:
    856         raise ValueError, "unknown algorithm '%s'"%algorithm
    857 
    858 def bessel_K(nu,z,algorithm="pari",prec=53):
    859     r"""
    860     Implements the "K-Bessel function", or "modified Bessel function,
    861     2nd kind", with index (or "order") nu and argument z. Defn::
    862    
    863                     pi*(bessel_I(-nu, z) - bessel_I(nu, z))
    864                    ----------------------------------------
    865                                 2*sin(pi*nu)
    866            
    867    
    868     if nu is not an integer and by taking a limit otherwise.
    869    
    870     Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
    871     PARI, nu can be complex and z must be real and positive.
    872    
    873     EXAMPLES::
    874    
    875         sage: bessel_K(3,2,"scipy")
    876         0.64738539094...
    877         sage: bessel_K(3,2)
    878         0.64738539094...
    879         sage: bessel_K(1,1)
    880         0.60190723019...
    881         sage: bessel_K(1,1,"pari",10)
    882         0.60
    883         sage: bessel_K(1,1,"pari",100)
    884         0.60190723019723457473754000154
    885 
    886     TESTS::
    887 
    888         sage: bessel_K(2,1.1, algorithm="maxima")
    889         Traceback (most recent call last):
    890         ...
    891         NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms
    892 
    893         Check whether the return value is real whenever the argument is real (#10251)::
    894 
    895         sage: bessel_K(5, 1.5, algorithm='scipy') in RR
    896         True
    897 
    898     """
    899     if algorithm=="scipy":
    900         if prec != 53:
    901             raise ValueError, "for the scipy algorithm the precision must be 53"
    902         import scipy.special
    903         ans = str(scipy.special.kv(float(nu),float(z)))
    904         ans = ans.replace("(","")
    905         ans = ans.replace(")","")
    906         ans = ans.replace("j","*I")
    907         ans = sage_eval(ans)
    908         return real(ans) if z in RR else ans
    909     elif algorithm == 'pari':
    910         from sage.libs.pari.all import pari
    911         try:
    912             R = RealField(prec)
    913             nu = R(nu)
    914             z = R(z)
    915         except TypeError:
    916             C = ComplexField(prec)
    917             nu = C(nu)
    918             z = C(z)
    919             K = C
    920         K = z.parent()
    921         return K(pari(nu).besselk(z, precision=prec))
    922     elif algorithm == 'maxima':
    923         raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
    924     else:
    925         raise ValueError, "unknown algorithm '%s'"%algorithm
    926    
    927 
    928 def bessel_Y(nu,z,algorithm="maxima", prec=53):
    929     r"""
    930     Implements the "Y-Bessel function", or "Bessel function of the 2nd
    931     kind", with index (or "order") nu and argument z.
    932    
    933     .. note::
    934 
    935        Currently only prec=53 is supported.
    936    
    937     Defn::
    938    
    939                     cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
    940                    -------------------------------------------------
    941                                      sin(nu*pi)
    942    
    943     if nu is not an integer and by taking a limit otherwise.
    944    
    945     Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
    946    
    947     This is computed using Maxima by default.
    948    
    949     EXAMPLES::
    950    
    951         sage: bessel_Y(2,1.1,"scipy")
    952         -1.4314714939...
    953         sage: bessel_Y(2,1.1)   
    954         -1.4314714939590...
    955         sage: bessel_Y(3.001,2.1)
    956         -1.0299574976424...
    957 
    958     TESTS::
    959 
    960         sage: bessel_Y(2,1.1, algorithm="pari")
    961         Traceback (most recent call last):
    962         ...
    963         NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
    964     """
    965     if algorithm=="scipy":
    966         if prec != 53:
    967             raise ValueError, "for the scipy algorithm the precision must be 53"
    968         import scipy.special
    969         ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
    970         ans = ans.replace("(","")
    971         ans = ans.replace(")","")
    972         ans = ans.replace("j","*I")
    973         ans = sage_eval(ans)
    974         return real(ans) if z in RR else ans
    975     elif algorithm == "maxima":
    976         if prec != 53:
    977             raise ValueError, "for the maxima algorithm the precision must be 53"
    978         return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    979     elif algorithm == "pari":
    980         raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
    981     else:
    982         raise ValueError, "unknown algorithm '%s'"%algorithm
    983    
    984 class Bessel():
    985     """
    986     A class implementing the I, J, K, and Y Bessel functions.
    987    
    988     EXAMPLES::
    989    
    990         sage: g = Bessel(2); g
    991         J_{2}
    992         sage: print g
    993         J-Bessel function of order 2
    994         sage: g.plot(0,10)
    995 
    996     ::
    997 
    998         sage: Bessel(2, typ='I')(pi)
    999         2.61849485263445
    1000         sage: Bessel(2, typ='J')(pi)
    1001         0.485433932631509
    1002         sage: Bessel(2, typ='K')(pi)
    1003         0.0510986902537926
    1004         sage: Bessel(2, typ='Y')(pi)
    1005         -0.0999007139289404
    1006     """
    1007     def __init__(self, nu, typ = "J", algorithm = None, prec = 53):
    1008         """
    1009         Initializes new instance of the Bessel class.
    1010 
    1011         INPUT:
    1012 
    1013          - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'
    1014            or 'Y'.
    1015 
    1016          - ``algorithm`` -- (default: maxima for type Y, pari for other types)
    1017            algorithm to use to compute the Bessel function: 'pari', 'maxima' or
    1018            'scipy'.  Note that type K is not implemented in Maxima and type Y
    1019            is not implemented in PARI.
    1020 
    1021          - ``prec`` -- (default: 53) precision in bits of the Bessel function.
    1022            Only supported for the PARI algorithm.
    1023 
    1024         EXAMPLES::
    1025 
    1026             sage: g = Bessel(2); g
    1027             J_{2}
    1028             sage: Bessel(1,'I')
    1029             I_{1}
    1030             sage: Bessel(6, prec=120)(pi)
    1031             0.014545966982505560573660369604001804
    1032             sage: Bessel(6, algorithm="pari")(pi)
    1033             0.0145459669825056
    1034 
    1035         For the Bessel J-function, Maxima returns a symbolic result.  For
    1036         types I and Y, we always get a numeric result::
    1037 
    1038             sage: b = Bessel(6, algorithm="maxima")(pi); b
    1039             bessel_j(6, pi)
    1040             sage: b.n(53)
    1041             0.0145459669825056
    1042             sage: Bessel(6, typ='I', algorithm="maxima")(pi)
    1043             0.0294619840059568
    1044             sage: Bessel(6, typ='Y', algorithm="maxima")(pi)
    1045             -4.33932818939038
    1046 
    1047         SciPy usually gives less precise results::
    1048 
    1049             sage: Bessel(6, algorithm="scipy")(pi)
    1050             0.0145459669825000...
    1051 
    1052         TESTS::
    1053 
    1054             sage: Bessel(1,'Z')
    1055             Traceback (most recent call last):
    1056             ...
    1057             ValueError: typ must be one of I, J, K, Y
    1058         """
    1059         if not (typ in ['I', 'J', 'K', 'Y']):
    1060             raise ValueError, "typ must be one of I, J, K, Y"
    1061 
    1062         # Did the user ask for the default algorithm?
    1063         if algorithm is None:
    1064             if typ == 'Y':
    1065                 algorithm = 'maxima'
    1066             else:
    1067                 algorithm = 'pari'
    1068 
    1069         self._system = algorithm
    1070         self._order = nu
    1071         self._type = typ
    1072         prec = int(prec)
    1073         if prec < 0:
    1074             raise ValueError, "prec must be a positive integer"
    1075         self._prec = int(prec)
    1076 
    1077     def __str__(self):
    1078         """
    1079         Returns a string representation of this Bessel object.
    1080 
    1081         TEST::
    1082 
    1083             sage: a = Bessel(1,'I')
    1084             sage: str(a)
    1085             'I-Bessel function of order 1'
    1086         """
    1087         return self.type()+"-Bessel function of order "+str(self.order())
    1088    
    1089     def __repr__(self):
    1090         """
    1091         Returns a string representation of this Bessel object.
    1092 
    1093         TESTS::
    1094 
    1095             sage: Bessel(1,'I')
    1096             I_{1}
    1097         """
    1098         return self.type()+"_{"+str(self.order())+"}"
    1099    
    1100     def type(self):
    1101         """
    1102         Returns the type of this Bessel object.
    1103 
    1104         TEST::
    1105 
    1106             sage: a = Bessel(3,'K')
    1107             sage: a.type()
    1108             'K'
    1109         """
    1110         return self._type
    1111    
    1112     def prec(self):
    1113         """
    1114         Returns the precision (in number of bits) used to represent this
    1115         Bessel function.
    1116 
    1117         TESTS::
    1118 
    1119             sage: a = Bessel(3,'K')
    1120             sage: a.prec()
    1121             53
    1122             sage: B = Bessel(20,prec=100); B
    1123             J_{20}
    1124             sage: B.prec()
    1125             100
    1126         """
    1127         return self._prec
    1128 
    1129     def order(self):
    1130         """
    1131         Returns the order of this Bessel function.
    1132 
    1133         TEST::
    1134 
    1135             sage: a = Bessel(3,'K')
    1136             sage: a.order()
    1137             3
    1138         """
    1139         return self._order
    1140 
    1141     def system(self):
    1142         """
    1143         Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with
    1144         this Bessel function.
    1145 
    1146         TESTS::
    1147 
    1148             sage: Bessel(20,algorithm='maxima').system()
    1149             'maxima'
    1150             sage: Bessel(20,prec=100).system()
    1151             'pari'
    1152         """
    1153         return self._system
    1154 
    1155     def __call__(self,z):
    1156         """
    1157         Implements evaluation of all the Bessel functions directly
    1158         from the Bessel class. This essentially allows a Bessel object to
    1159         behave like a function that can be invoked.
    1160 
    1161         TESTS::
    1162 
    1163             sage: Bessel(3,'K')(5.0)
    1164             0.00829176841523093
    1165             sage: Bessel(20,algorithm='maxima')(5.0)
    1166             2.77033005213e-11
    1167             sage: Bessel(20,prec=100)(5.0101010101010101)
    1168             2.8809188227195382093062257967e-11
    1169             sage: B = Bessel(2,'Y',algorithm='scipy',prec=50)
    1170             sage: B(2.0)
    1171             Traceback (most recent call last):
    1172             ...
    1173             ValueError: for the scipy algorithm the precision must be 53
    1174         """
    1175         nu = self.order()
    1176         t = self.type()
    1177         s = self.system()
    1178         p = self.prec()
    1179         if t == "I":
    1180             return bessel_I(nu,z,algorithm=s,prec=p)
    1181         if t == "J":
    1182             return bessel_J(nu,z,algorithm=s,prec=p)
    1183         if t == "K":
    1184             return bessel_K(nu,z,algorithm=s,prec=p)
    1185         if t == "Y":
    1186             return bessel_Y(nu,z,algorithm=s,prec=p)
    1187        
    1188     def plot(self,a,b):
    1189         """
    1190         Enables easy plotting of all the Bessel functions directly
    1191         from the Bessel class.
    1192 
    1193         TESTS::
    1194 
    1195             sage: plot(Bessel(2),3,4)
    1196             sage: Bessel(2).plot(3,4)
    1197             sage: P = Bessel(2,'I').plot(1,5)
    1198             sage: P += Bessel(2,'J').plot(1,5)
    1199             sage: P += Bessel(2,'K').plot(1,5)
    1200             sage: P += Bessel(2,'Y').plot(1,5)
    1201             sage: show(P)
    1202         """
    1203         nu = self.order()
    1204         s = self.system()
    1205         t = self.type()
    1206         if t == "I":
    1207             f = lambda z: bessel_I(nu,z,s) 
    1208             P = plot(f,a,b)
    1209         if t == "J":
    1210             f = lambda z: bessel_J(nu,z,s)
    1211             P = plot(f,a,b)
    1212         if t == "K":
    1213             f = lambda z: bessel_K(nu,z,s)
    1214             P = plot(f,a,b)
    1215         if t == "Y":
    1216             f = lambda z: bessel_Y(nu,z,s)
    1217             P = plot(f,a,b)
    1218         return P
    1219    
    1220620def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53):
    1221621    r"""
    1222622    Default is a wrap of PARI's hyperu(alpha,beta,x) function.