Ticket #4102: trac_symbolic_bessel_v2.patch

File trac_symbolic_bessel_v2.patch, 43.6 KB (added by benjaminfjones, 8 years ago)

more work in progress

  • sage/functions/all.py

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1356654134 28800
    # Node ID e893ca8c441435d3adad6d3c8914af15680ce560
    # Parent  9519a7bb2f424429ef98c756b2f4d0ea2fec8c8f
    Trac 4102: Implement symbolic Bessel functions
    
    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 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"""
     2This module provides symbolic Bessel Functions.
     3
     4AUTHORS:
     5
     6    - Benjamin Jones (2012-12-27): initial version
     7
     8This module provides symbolic Bessel functions.
     9
     10-  Bessel functions, first defined by the Swiss mathematician
     11   Daniel Bernoulli and named after Friedrich Bessel, are canonical
     12   solutions y(x) of Bessel's differential equation:
     13
     14
     15   .. math::
     16
     17         x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0,
     18
     19
     20   for an arbitrary real number `\alpha` (the order).
     21
     22-  Another important formulation of the two linearly independent
     23   solutions to Bessel's equation are the Hankel functions
     24   `H_\alpha^{(1)}(x)` and `H_\alpha^{(2)}(x)`,
     25   defined by:
     26
     27
     28   .. math::
     29
     30         H_\alpha^{(1)}(x) = J_\alpha(x) + i Y_\alpha(x)
     31
     32
     33
     34   .. math::
     35
     36         H_\alpha^{(2)}(x) = J_\alpha(x) - i Y_\alpha(x)
     37
     38   where `i` is the imaginary unit (and `J_*` and
     39   `Y_*` are the usual J- and Y-Bessel functions). These
     40   linear combinations are also known as Bessel functions of the third
     41   kind; they are two linearly independent solutions of Bessel's
     42   differential equation. They are named for Hermann Hankel.
     43
     44REFERENCES:
     45
     46    - Abramowitz and Stegun: Handbook of Mathematical Functions,
     47      http://www.math.sfu.ca/~cbm/aands/
     48
     49    - http://en.wikipedia.org/wiki/Bessel_function
     50"""
     51
     52#*****************************************************************************
     53#       Copyright (C) 2012 Benjamin Jones <benjaminfjones@gmail.com>
     54#
     55#  Distributed under the terms of the GNU General Public License (GPL)
     56#
     57#    This code is distributed in the hope that it will be useful,
     58#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     59#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     60#    General Public License for more details.
     61#
     62#  The full text of the GPL is available at:
     63#
     64#                  http://www.gnu.org/licenses/
     65#*****************************************************************************
     66
     67from sage.rings.all import RR
     68from sage.symbolic.function import BuiltinFunction, is_inexact
     69from sage.symbolic.expression import Expression
     70from sage.structure.coerce import parent
     71import sage.structure.element
     72from sage.calculus.calculus import maxima
     73from sage.libs.mpmath import utils as mpmath_utils
     74
     75class Function_Bessel_I(BuiltinFunction):
     76    r"""
     77    The Bessel I function.
     78
     79    DEFINITION:
     80
     81    .. math:
     82
     83        bessel\_I(n, z) = I_n(z) = ...
     84
     85    The derivative with respect to `z` is:
     86
     87    ..math::
     88
     89        \frac{d}{dz} I_n(z) = \frac{1}{z^n} \left\(z^n I_{n-1}(z) - n z^{n-1} I_n(z) \right\)
     90
     91    EXAMPLES::
     92
     93        sage: bessel_I(1, x)
     94        bessel_I(1, x)
     95        sage: bessel_I(1.0, 1.0)
     96        0.565159103992485
     97        sage: n = var('n')
     98        sage: bessel_I(n, x)
     99        bessel_I(n, x)
     100        sage: bessel_I(2, I).n(digits=30)
     101        -0.114903484931900480469646881335
     102
     103    Examples of symbolic manipulation::
     104
     105        sage: a = bessel_I(pi, bessel_I(1, I))
     106        sage: N(a, digits=20)
     107        0.00026073272117205890528 - 0.0011528954889080572266*I
     108
     109        sage: f = bessel_I(2, x)
     110        sage: f.diff(x)
     111        -2*bessel_I(2, x)/x + bessel_I(1, x)
     112
     113    Visualization::
     114
     115        sage: plot(bessel_I(1,x), (x,0,5), color='blue')
     116        sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5)) # long time
     117
     118    ALGORITHM:
     119
     120        Numerical evaluation is handled by the mpmath library. Symbolics are
     121        handled by a combination of Maxima and Sage (Ginac/Pynac).
     122    """
     123    def __init__(self):
     124        """
     125        See the docstring for :meth:`Function_Bessel_I`.
     126
     127        EXAMPLES::
     128
     129            sage: bessel_I(1,x)
     130            bessel_I(1, x)
     131        """
     132        BuiltinFunction.__init__(self, "bessel_I", nargs=2,
     133                                 latex_name=r'bessel_I',
     134                                 conversions=dict(maxima='bessel_i', mathematica='BesselI'))
     135
     136    def _eval_(self, n, z):
     137        """
     138        EXAMPLES::
     139
     140            sage: y=var('y')
     141            sage: bessel_I(y,x)
     142            bessel_I(y, x)
     143            sage: bessel_I(0.0, 1.0)
     144            1.26606587775201
     145        """
     146        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
     147               (is_inexact(n) or is_inexact(z)):
     148            coercion_model = sage.structure.element.get_coercion_model()
     149            n, z = coercion_model.canonical_coercion(n, z)
     150            return self._evalf_(n, z, parent(n))
     151
     152        return None # leaves the expression unevaluated
     153
     154    def _evalf_(self, n, z, parent=None):
     155        """
     156        EXAMPLES::
     157
     158            sage: bessel_I(1,3).n(digits=20)
     159            3.9533702174026093965
     160        """
     161        import mpmath
     162        return mpmath_utils.call(mpmath.besseli, n, z, parent=parent)
     163
     164    def _derivative_(self, n, z, diff_param=None):
     165        """
     166        Return the derivative of the Bessel I function.
     167
     168        EXAMPLES::
     169
     170            sage: f(z) = bessel_I(10, z)
     171            sage: derivative(f, z)
     172            z |--> -10*bessel_I(10, z)/z + bessel_I(9, z)
     173        """
     174        return bessel_I(n-1, z) - n/z*bessel_I(n, z)
     175
     176bessel_I = Function_Bessel_I()
     177
     178class Function_Bessel_J(BuiltinFunction):
     179    r"""
     180    The Bessel J function.
     181
     182    The derivative with respect to `z` is:
     183
     184    ..math::
     185
     186        \frac{d}{dz} J_n(z) = \frac{1}{z^n} \left\(z^n J_{n-1}(z) - n z^{n-1} J_n(z) \right\)
     187
     188    EXAMPLES::
     189
     190        sage: bessel_J(1, x)
     191        bessel_J(1, x)
     192        sage: bessel_J(1.0, 1.0)
     193        0.440050585744933
     194        sage: n = var('n')
     195        sage: bessel_J(n, x)
     196        bessel_J(n, x)
     197        sage: bessel_J(2, I).n(digits=30)
     198        -0.135747669767038281182852569995
     199
     200    Examples of symbolic manipulation::
     201
     202        sage: a = bessel_J(pi, bessel_J(1, I))
     203        sage: N(a, digits=20)
     204        0.00059023706363796717363 - 0.0026098820470081958110*I
     205
     206        sage: f = bessel_J(2, x)
     207        sage: f.diff(x)
     208        -2*bessel_J(2, x)/x + bessel_J(1, x)
     209
     210    Currently, integration is not supported (directly) since we cannot yet convert
     211    hypergeometric functions to and from Maxima::
     212
     213        sage: f = bessel_J(2, x)
     214        sage: f.integrate(x)
     215        Traceback (most recent call last):
     216        ...
     217        TypeError: cannot coerce arguments: no canonical coercion from <type 'list'> to Symbolic Ring
     218
     219        sage: m = maxima(bessel_J(2, x))
     220        sage: m.integrate(x)
     221        hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24
     222
     223    Visualization::
     224
     225        sage: plot(bessel_J(1,x), (x,0,5), color='blue')
     226        sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5)) # long time
     227
     228    ALGORITHM:
     229
     230        Numerical evaluation is handled by the mpmath library. Symbolics are
     231        handled by a combination of Maxima and Sage (Ginac/Pynac).
     232    """
     233    def __init__(self):
     234        """
     235        See the docstring for :meth:`Function_Bessel_J`.
     236        """
     237        BuiltinFunction.__init__(self, "bessel_J", nargs=2,
     238                                 latex_name=r'bessel_J',
     239                                 conversions=dict(maxima='bessel_j', mathematica='BesselJ'))
     240
     241    def _eval_(self, n, z):
     242        """
     243        """
     244        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
     245               (is_inexact(n) or is_inexact(z)):
     246            coercion_model = sage.structure.element.get_coercion_model()
     247            n, z = coercion_model.canonical_coercion(n, z)
     248            return self._evalf_(n, z, parent(n))
     249
     250        return None # leaves the expression unevaluated
     251
     252    def _evalf_(self, n, z, parent=None):
     253        """
     254        """
     255        import mpmath
     256        return mpmath_utils.call(mpmath.besselj, n, z, parent=parent)
     257
     258    def _derivative_(self, n, z, diff_param=None):
     259        """
     260        Return the derivative of the Bessel I function.
     261
     262        EXAMPLES::
     263
     264            sage: f(z) = bessel_J(10, z)
     265            sage: derivative(f, z)
     266            z |--> -10*bessel_J(10, z)/z + bessel_J(9, z)
     267        """
     268        return bessel_J(n-1, z) - n/z*bessel_J(n, z)
     269
     270bessel_J = Function_Bessel_J()
     271
     272class Function_Bessel_K(BuiltinFunction):
     273    r"""
     274    The Bessel K function.
     275
     276    DEFINITION:
     277
     278    .. math:
     279
     280        bessel\_K(n, z) = K_n(z) = ...
     281
     282    The derivative with respect to `z` is:
     283
     284    ..math::
     285
     286        \frac{d}{dz} K_n(z) = \frac{1}{z^n} \left\(z^n K_{n-1}(z) - n z^{n-1} K_n(z) \right\)
     287
     288    """
     289    def __init__(self):
     290        """
     291        See the docstring for :meth:`Function_Bessel_K`.
     292        """
     293        BuiltinFunction.__init__(self, "bessel_K", nargs=2,
     294                                 latex_name=r'bessel_K',
     295                                 conversions=dict(maxima='bessel_k', mathematica='BesselK'))
     296
     297    def _eval_(self, n, z):
     298        """
     299        """
     300        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
     301               (is_inexact(n) or is_inexact(z)):
     302            coercion_model = sage.structure.element.get_coercion_model()
     303            n, z = coercion_model.canonical_coercion(n, z)
     304            return self._evalf_(n, z, parent(n))
     305
     306        return None # leaves the expression unevaluated
     307
     308    def _evalf_(self, n, z, parent=None):
     309        """
     310        """
     311        import mpmath
     312        return mpmath_utils.call(mpmath.besselk, n, z, parent=parent)
     313
     314    def _derivative_(self, n, z, diff_param=None):
     315        """
     316        Return the derivative of the Bessel K function.
     317
     318        EXAMPLES::
     319
     320            sage: f(z) = bessel_K(10, z)
     321            sage: derivative(f, z)
     322            z |--> -10*bessel_K(10, z)/z + bessel_K(9, z)
     323        """
     324        return bessel_K(n-1, z) - n/z*bessel_K(n, z)
     325
     326bessel_K = Function_Bessel_K()
     327
     328class Function_Bessel_Y(BuiltinFunction):
     329    r"""
     330    The Bessel Y function.
     331
     332    DEFINITION:
     333
     334    .. math:
     335
     336        bessel\_Y(n, z) = Y_n(z) = ...
     337
     338    The derivative with respect to `z` is:
     339
     340    ..math::
     341
     342        \frac{d}{dz} Y_n(z) = \frac{1}{z^n} \left\(z^n Y_{n-1}(z) - n z^{n-1} Y_n(z) \right\)
     343
     344    """
     345    def __init__(self):
     346        """
     347        See the docstring for :meth:`Function_Bessel_Y`.
     348        """
     349        BuiltinFunction.__init__(self, "bessel_Y", nargs=2,
     350                                 latex_name=r'bessel_Y',
     351                                 conversions=dict(maxima='bessel_y', mathematica='BesselY'))
     352
     353    def _eval_(self, n, z):
     354        """
     355        """
     356        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
     357               (is_inexact(n) or is_inexact(z)):
     358            coercion_model = sage.structure.element.get_coercion_model()
     359            n, z = coercion_model.canonical_coercion(n, z)
     360            return self._evalf_(n, z, parent(n))
     361
     362        return None # leaves the expression unevaluated
     363
     364    def _evalf_(self, n, z, parent=None):
     365        """
     366        """
     367        import mpmath
     368        return mpmath_utils.call(mpmath.bessely, n, z, parent=parent)
     369
     370    def _derivative_(self, n, z, diff_param=None):
     371        """
     372        Return the derivative of the Bessel Y function.
     373
     374        EXAMPLES::
     375
     376            sage: f(z) = bessel_Y(10, z)
     377            sage: derivative(f, z)
     378            z |--> -10*bessel_Y(10, z)/z + bessel_Y(9, z)
     379        """
     380        return bessel_Y(n-1, z) - n/z*bessel_Y(n, z)
     381
     382bessel_Y = Function_Bessel_Y()
     383
     384# dictionary used in Bessel
     385bessel_type_dict = {'I': bessel_I, 'J': bessel_J, 'K': bessel_K, 'Y': bessel_Y}
     386
     387def Bessel(*args, **kwds):
     388    """
     389    A factory function that produces symbolic I, J, K, and Y Bessel functions. There
     390    are several ways to call this function:
     391
     392        - ``Bessel(order, type)``
     393        - ``Bessel(order)`` -- type defaults to 'J'
     394        - ``Bessel(order, typ=T)``
     395        - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter function
     396        - ``Bessel()`` -- order is unspecified, type is 'J'
     397
     398    where `order` can be any integer and T must be one of the strings 'I', 'J',
     399    'K', or 'Y'.
     400
     401    See the EXAMPLES below.
     402
     403    EXAMPLES:
     404
     405    Construction of Bessel functions with various orders and types::
     406
     407        sage: Bessel()
     408        bessel_J
     409        sage: Bessel(1)(x)
     410        bessel_J(1, x)
     411        sage: Bessel(1, 'Y')(x)
     412        bessel_Y(1, x)
     413        sage: Bessel(-2, 'Y')(x)
     414        bessel_Y(-2, x)
     415        sage: Bessel(typ='K')
     416        bessel_K
     417        sage: Bessel(0, typ='I')(x)
     418        bessel_I(0, x)
     419
     420    Evaluation::
     421
     422        sage: f = Bessel(1)
     423        sage: f(3.0)
     424        0.339058958525936
     425        sage: f(3)
     426        bessel_J(1, 3)
     427        sage: f(3).n(digits=50)
     428        0.33905895852593645892551459720647889697308041819801
     429
     430        sage: g = Bessel(typ='J')
     431        sage: g(1,3)
     432        bessel_J(1, 3)
     433        sage: g(2, 3+I).n()
     434        0.634160370148554 + 0.0253384000032695*I
     435        sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15
     436        True
     437
     438    Symbolic calculus::
     439
     440        sage: f(x) = Bessel(0, 'J')(x)
     441        sage: derivative(f, x)
     442        x |--> bessel_J(-1, x)
     443        sage: derivative(f, x, x)
     444        x |--> bessel_J(-1, x)/x + bessel_J(-2, x)
     445
     446    Conversion to other systems::
     447
     448        sage: x,y = var('x,y')
     449        sage: f = maxima(Bessel(typ='K')(x,y))
     450        sage: f.derivative('x')
     451        %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)
     452        sage: f.derivative('y')
     453        -(bessel_k(x+1,y)+bessel_k(x-1,y))/2
     454
     455    Plotting::
     456
     457        sage: f(x) = Bessel(0)(x); f
     458        x |--> bessel_J(0, x)
     459        sage: plot(f, (x, 1, 10))
     460
     461        sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10)
     462
     463        sage: G = Graphics()
     464        sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ])
     465        sage: show(G)
     466
     467    A recreation of Abramowitz and Stegun Figure 9.1::
     468
     469        sage: G  = plot(Bessel(0, 'J'), 0, 15, color='black')
     470        sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black')
     471        sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted')
     472        sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted')
     473        sage: show(G, ymin=-1, ymax=1)
     474
     475    """
     476    # Determine the order and type of function from the arguments and keywords.
     477    # These are recored in local variables: _type, _order, _system, _nargs.
     478    _type = None
     479    if len(args) == 0: # no order specified
     480        _order = None
     481        _nargs = 2
     482    elif len(args) == 1: # order is specified
     483        _order = args[0]
     484        _nargs = 1
     485    elif len(args) == 2: # both order and type are positional arguments
     486        _order = args[0]
     487        _type = args[1]
     488        _nargs = 1
     489    else:
     490        raise TypeError("at most two positional arguments may be specified, "
     491                       +"see the docstring for Bessel")
     492    # check for type inconsistency
     493    if _type is not None and 'typ' in kwds and _type != kwds['typ']:
     494        raise ValueError("inconsistent types given")
     495    # record the function type
     496    if _type is None:
     497        if 'typ' in kwds:
     498            _type = kwds['typ']
     499        else:
     500            _type = 'J'
     501    if not (_type in ['I', 'J', 'K', 'Y']):
     502        raise ValueError("type must be one of I, J, K, Y")
     503    # record the numerical evaluation system
     504    if 'algorithm' in kwds:
     505        _system = kwds['algorithm']
     506    else:
     507        _system = 'mpmath'
     508
     509    # return the function
     510    # TODO remove assertions
     511    assert _type in ['I', 'J', 'K', 'Y']
     512    assert _order is None or _order in RR
     513    assert _nargs == 1 or _nargs == 2
     514    assert _system == 'mpmath'
     515
     516    # TODO what to do with _system?
     517    _f = bessel_type_dict[_type]
     518    if _nargs == 1:
     519        return lambda x: _f(_order, x)
     520    else:
     521        return _f
  • 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
    39 
    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 
    7532-  Airy function The function `Ai(x)` and the related
    7633   function `Bi(x)`, which is also called an Airy function,
    7734   are solutions to the differential equation
    7835
    79    
     36
    8037   .. math::
    8138
    82          y'' - xy = 0, 
     39         y'' - xy = 0,
    8340
    8441   known as the Airy equation. They belong to the class of 'Bessel functions of
    8542   fractional order'. The initial conditions
     
    9350
    9451-  Spherical harmonics: Laplace's equation in spherical coordinates
    9552   is:
    96    
     53
    9754   .. math::
    9855
    99        {\frac{1}{r^2}}{\frac{\partial}{\partial r}}   \left(r^2 {\frac{\partial f}{\partial r}}\right) +   {\frac{1}{r^2}\sin\theta}{\frac{\partial}{\partial \theta}}   \left(\sin\theta {\frac{\partial f}{\partial \theta}}\right) +   {\frac{1}{r^2\sin^2\theta}}{\frac{\partial^2 f}{\partial \varphi^2}} = 0. 
     56       {\frac{1}{r^2}}{\frac{\partial}{\partial r}}   \left(r^2 {\frac{\partial f}{\partial r}}\right) +   {\frac{1}{r^2}\sin\theta}{\frac{\partial}{\partial \theta}}   \left(\sin\theta {\frac{\partial f}{\partial \theta}}\right) +   {\frac{1}{r^2\sin^2\theta}}{\frac{\partial^2 f}{\partial \varphi^2}} = 0.
    10057
    10158
    10259   Note that the spherical coordinates `\theta` and
     
    10764
    10865   The general solution which remains finite towards infinity is a
    10966   linear combination of functions of the form
    110    
     67
    11168   .. math::
    11269
    113          r^{-1-\ell} \cos (m \varphi) P_\ell^m (\cos{\theta} ) 
     70         r^{-1-\ell} \cos (m \varphi) P_\ell^m (\cos{\theta} )
    11471
    11572
    11673   and
    117    
     74
    11875   .. math::
    11976
    120          r^{-1-\ell} \sin (m \varphi) P_\ell^m (\cos{\theta} ) 
     77         r^{-1-\ell} \sin (m \varphi) P_\ell^m (\cos{\theta} )
    12178
    12279
    12380   where `P_\ell^m` are the associated Legendre polynomials,
     
    12683   solutions with integer parameters `\ell \ge 0` and
    12784   `- \ell\leq m\leq \ell`, can be written as linear
    12885   combinations of:
    129    
     86
    13087   .. math::
    13188
    132          U_{\ell,m}(r,\theta , \varphi ) = r^{-1-\ell} Y_\ell^m( \theta , \varphi ) 
     89         U_{\ell,m}(r,\theta , \varphi ) = r^{-1-\ell} Y_\ell^m( \theta , \varphi )
    13390
    13491
    13592   where the functions `Y` are the spherical harmonic
    13693   functions with parameters `\ell`, `m`, which can be
    13794   written as:
    138    
     95
    13996   .. math::
    14097
    141          Y_\ell^m( \theta , \varphi )     = \sqrt{{\frac{(2\ell+1)}{4\pi}}{\frac{(\ell-m)!}{(\ell+m)!}}}       \cdot e^{i m \varphi } \cdot P_\ell^m ( \cos{\theta} ) . 
     98         Y_\ell^m( \theta , \varphi )     = \sqrt{{\frac{(2\ell+1)}{4\pi}}{\frac{(\ell-m)!}{(\ell+m)!}}}       \cdot e^{i m \varphi } \cdot P_\ell^m ( \cos{\theta} ) .
    14299
    143100
    144101
    145102   The spherical harmonics obey the normalisation condition
    146103
    147    
     104
    148105   .. math::
    149106
    150      \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta . 
     107     \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta .
    151108
    152109
    153110
    154111-  When solving for separable solutions of Laplace's equation in
    155112   spherical coordinates, the radial equation has the form:
    156    
     113
    157114   .. math::
    158115
    159          x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0. 
     116         x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0.
    160117
    161118
    162119   The spherical Bessel functions `j_n` and `y_n`,
    163120   are two linearly independent solutions to this equation. They are
    164121   related to the ordinary Bessel functions `J_n` and
    165122   `Y_n` by:
    166    
     123
    167124   .. math::
    168125
    169          j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x), 
     126         j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x),
    170127
    171128
    172    
     129
    173130   .. math::
    174131
    175          y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x)     = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x). 
     132         y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x)     = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x).
    176133
    177134
    178135
     
    180137   `y = U(a,b,x)` is defined to be the solution to Kummer's
    181138   differential equation
    182139
    183    
     140
    184141   .. math::
    185142
    186      xy'' + (b-x)y' - ay = 0, 
     143     xy'' + (b-x)y' - ay = 0,
    187144
    188145   which satisfies `U(a,b,x) \sim x^{-a}`, as
    189146   `x\rightarrow \infty`. (There is a linearly independent
     
    213170   doubly-periodic, meromorphic functions satisfying the following
    214171   three properties:
    215172
    216    
     173
    217174   #. There is a simple zero at the corner p, and a simple pole at the
    218175      corner q.
    219176
     
    232189      is 1.
    233190
    234191
    235    We can write 
    236    
     192   We can write
     193
    237194   .. math::
    238195
    239       pq(x)=\frac{pr(x)}{qr(x)} 
     196      pq(x)=\frac{pr(x)}{qr(x)}
    240197
    241198
    242199   where `p`, `q`, and `r` are any of the
     
    244201   the understanding that `ss=cc=dd=nn=1`.
    245202
    246203   Let
    247    
     204
    248205   .. math::
    249206
    250          u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}} 
     207         u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}}
    251208
    252209
    253210   Then the *Jacobi elliptic function* `sn(u)` is given by
    254    
     211
    255212   .. math::
    256213
    257         {sn}\; u = \sin \phi 
     214        {sn}\; u = \sin \phi
    258215
    259216   and `cn(u)` is given by
    260    
     217
    261218   .. math::
    262219
    263        {cn}\; u = \cos \phi 
     220       {cn}\; u = \cos \phi
    264221
    265222   and
    266    
     223
    267224   .. math::
    268225
    269      {dn}\; u = \sqrt {1-m\sin^2 \phi}. 
     226     {dn}\; u = \sqrt {1-m\sin^2 \phi}.
    270227
    271228   To emphasize the dependence on `m`, one can write
    272229   `sn(u,m)` for example (and similarly for `cn` and
     
    276233   solutions to the following nonlinear ordinary differential
    277234   equations:
    278235
    279    
     236
    280237   -  `\mathrm{sn}\,(x;k)` solves the differential equations
    281      
     238
    282239      .. math::
    283240
    284241          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1+k^2) y - 2 k^2 y^3 = 0,
    285          
     242
    286243      and
    287244
    288245      .. math::
     
    291248
    292249   -  `\mathrm{cn}\,(x;k)` solves the differential equations
    293250
    294      
     251
    295252      .. math::
    296      
    297          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0, 
     253
     254         \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0,
    298255
    299256      and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 + k^2 y^2)`.
    300257
    301258   -  `\mathrm{dn}\,(x;k)` solves the differential equations
    302      
     259
    303260      .. math::
    304261
    305          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0, 
     262         \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0,
    306263
    307264      and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2= y^2 (1 - k^2 - y^2)`.
    308265
     
    320277      `dn(x, 1) = \mbox{\rm sech} x`.
    321278
    322279   -  The incomplete elliptic integrals (of the first kind, etc.) are:
    323      
     280
    324281      .. math::
    325  
    326          \begin{array}{c} \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2}}\, dx,\\ \displaystyle\int_0^\phi \sqrt{1 - m\sin(x)^2}\, dx,\\ \displaystyle\int_0^\phi \frac{\sqrt{1-mt^2}}{\sqrt(1 - t^2)}\, dx,\\ \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2\sqrt{1 - n\sin(x)^2}}}\, dx, \end{array} 
     282
     283         \begin{array}{c} \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2}}\, dx,\\ \displaystyle\int_0^\phi \sqrt{1 - m\sin(x)^2}\, dx,\\ \displaystyle\int_0^\phi \frac{\sqrt{1-mt^2}}{\sqrt(1 - t^2)}\, dx,\\ \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2\sqrt{1 - n\sin(x)^2}}}\, dx, \end{array}
    327284
    328285      and the complete ones are obtained by taking `\phi =\pi/2`.
    329286
    330        
    331287REFERENCES:
    332288
    333289- Abramowitz and Stegun: Handbook of Mathematical Functions,
    334290  http://www.math.sfu.ca/~cbm/aands/
    335291
    336 - http://en.wikipedia.org/wiki/Bessel_function
    337 
    338292- http://en.wikipedia.org/wiki/Airy_function
    339293
    340294- http://en.wikipedia.org/wiki/Spherical_harmonics
     
    349303- Online Encyclopedia of Special Function
    350304  http://algo.inria.fr/esf/index.html
    351305
    352 TODO: Resolve weird bug in commented out code in hypergeometric_U
    353 below.
    354 
    355306AUTHORS:
    356307
    357308- David Joyner and William Stein
     
    652603   _init()
    653604   return RDF(meval("airy_bi(%s)"%RDF(x)))
    654605
    655 
    656 def bessel_I(nu,z,algorithm = "pari",prec=53):
    657     r"""
    658     Implements the "I-Bessel function", or "modified Bessel function,
    659     1st kind", with index (or "order") nu and argument z.
    660    
    661     INPUT:
    662    
    663    
    664     -  ``nu`` - a real (or complex, for pari) number
    665    
    666     -  ``z`` - a real (positive) algorithm - "pari" or
    667        "maxima" or "scipy" prec - real precision (for PARI only)
    668    
    669    
    670     DEFINITION::
    671    
    672             Maxima:
    673                              inf
    674                             ====   - nu - 2 k  nu + 2 k
    675                             \     2          z
    676                              >    -------------------
    677                             /     k! Gamma(nu + k + 1)
    678                             ====
    679                             k = 0
    680        
    681             PARI:
    682            
    683                              inf
    684                             ====   - 2 k  2 k
    685                             \     2      z    Gamma(nu + 1)
    686                              >    -----------------------
    687                             /       k! Gamma(nu + k + 1)
    688                             ====
    689                             k = 0
    690        
    691            
    692    
    693     Sometimes ``bessel_I(nu,z)`` is denoted
    694     ``I_nu(z)`` in the literature.
    695    
    696     .. warning::
    697 
    698        In Maxima (the manual says) i0 is deprecated but
    699        ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
    700        though.)
    701    
    702     EXAMPLES::
    703    
    704         sage: bessel_I(1,1,"pari",500)
    705         0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
    706         sage: bessel_I(1,1)
    707         0.565159103992485
    708         sage: bessel_I(2,1.1,"maxima") 
    709         0.16708949925104...
    710         sage: bessel_I(0,1.1,"maxima")
    711         1.32616018371265...
    712         sage: bessel_I(0,1,"maxima")   
    713         1.2660658777520...
    714         sage: bessel_I(1,1,"scipy")
    715         0.565159103992...
    716 
    717     Check whether the return value is real whenever the argument is real (#10251)::
    718    
    719         sage: bessel_I(5, 1.5, algorithm='scipy') in RR
    720         True
    721        
    722     """
    723     if algorithm=="pari":
    724         from sage.libs.pari.all import pari
    725         try:
    726             R = RealField(prec)
    727             nu = R(nu)
    728             z = R(z)
    729         except TypeError:
    730             C = ComplexField(prec)
    731             nu = C(nu)
    732             z = C(z)
    733             K = C
    734         K = z.parent()
    735         return K(pari(nu).besseli(z, precision=prec))
    736     elif algorithm=="scipy":
    737         if prec != 53:
    738             raise ValueError, "for the scipy algorithm the precision must be 53"
    739         import scipy.special
    740         ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
    741         ans = ans.replace("(","")
    742         ans = ans.replace(")","")
    743         ans = ans.replace("j","*I")
    744         ans = sage_eval(ans)
    745         return real(ans) if z in RR else ans # Return real value when arg is real
    746     elif algorithm == "maxima":
    747         if prec != 53:
    748             raise ValueError, "for the maxima algorithm the precision must be 53"
    749         return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
    750     else:
    751         raise ValueError, "unknown algorithm '%s'"%algorithm
    752        
    753 def bessel_J(nu,z,algorithm="pari",prec=53):
    754     r"""
    755     Return value of the "J-Bessel function", or "Bessel function, 1st
    756     kind", with index (or "order") nu and argument z.
    757    
    758     ::
    759    
    760             Defn:
    761             Maxima:
    762                              inf
    763                             ====          - nu - 2 k  nu + 2 k
    764                             \     (-1)^k 2           z
    765                              >    -------------------------
    766                             /        k! Gamma(nu + k + 1)
    767                             ====
    768                             k = 0
    769        
    770             PARI:
    771            
    772                              inf
    773                             ====          - 2k    2k
    774                             \     (-1)^k 2      z    Gamma(nu + 1)
    775                              >    ----------------------------
    776                             /         k! Gamma(nu + k + 1)
    777                             ====
    778                             k = 0
    779            
    780    
    781     Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
    782    
    783     .. warning::
    784 
    785        Inaccurate for small values of z.
    786    
    787     EXAMPLES::
    788    
    789         sage: bessel_J(2,1.1)
    790         0.136564153956658
    791         sage: bessel_J(0,1.1)
    792         0.719622018527511
    793         sage: bessel_J(0,1)
    794         0.765197686557967
    795         sage: bessel_J(0,0)
    796         1.00000000000000
    797         sage: bessel_J(0.1,0.1)
    798         0.777264368097005
    799    
    800     We check consistency of PARI and Maxima::
    801    
    802         sage: n(bessel_J(3,10,"maxima"))
    803         0.0583793793051...
    804         sage: n(bessel_J(3,10,"pari")) 
    805         0.0583793793051868
    806         sage: bessel_J(3,10,"scipy")
    807         0.0583793793052...
    808 
    809     Check whether the return value is real whenever the argument is real (#10251)::                                                                                                                                                           
    810         sage: bessel_J(5, 1.5, algorithm='scipy') in RR                                                                     
    811         True
    812     """
    813    
    814     if algorithm=="pari":
    815         from sage.libs.pari.all import pari
    816         try:
    817             R = RealField(prec)
    818             nu = R(nu)
    819             z = R(z)
    820         except TypeError:
    821             C = ComplexField(prec)
    822             nu = C(nu)
    823             z = C(z)
    824             K = C
    825         if nu == 0:
    826             nu = ZZ(0)
    827         K = z.parent()
    828         return K(pari(nu).besselj(z, precision=prec))
    829     elif algorithm=="scipy":
    830         if prec != 53:
    831             raise ValueError, "for the scipy algorithm the precision must be 53"
    832         import scipy.special
    833         ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
    834         ans = ans.replace("(","")
    835         ans = ans.replace(")","")
    836         ans = ans.replace("j","*I")
    837         ans = sage_eval(ans)
    838         return real(ans) if z in RR else ans
    839     elif algorithm == "maxima":
    840         if prec != 53:
    841             raise ValueError, "for the maxima algorithm the precision must be 53"
    842         return maxima_function("bessel_j")(nu, z)
    843     else:
    844         raise ValueError, "unknown algorithm '%s'"%algorithm
    845 
    846 def bessel_K(nu,z,algorithm="pari",prec=53):
    847     r"""
    848     Implements the "K-Bessel function", or "modified Bessel function,
    849     2nd kind", with index (or "order") nu and argument z. Defn::
    850    
    851                     pi*(bessel_I(-nu, z) - bessel_I(nu, z))
    852                    ----------------------------------------
    853                                 2*sin(pi*nu)
    854            
    855    
    856     if nu is not an integer and by taking a limit otherwise.
    857    
    858     Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
    859     PARI, nu can be complex and z must be real and positive.
    860    
    861     EXAMPLES::
    862    
    863         sage: bessel_K(3,2,"scipy")
    864         0.64738539094...
    865         sage: bessel_K(3,2)
    866         0.64738539094...
    867         sage: bessel_K(1,1)
    868         0.60190723019...
    869         sage: bessel_K(1,1,"pari",10)
    870         0.60
    871         sage: bessel_K(1,1,"pari",100)
    872         0.60190723019723457473754000154
    873 
    874     TESTS::
    875 
    876         sage: bessel_K(2,1.1, algorithm="maxima")
    877         Traceback (most recent call last):
    878         ...
    879         NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms
    880 
    881         Check whether the return value is real whenever the argument is real (#10251)::
    882 
    883         sage: bessel_K(5, 1.5, algorithm='scipy') in RR
    884         True
    885 
    886     """
    887     if algorithm=="scipy":
    888         if prec != 53:
    889             raise ValueError, "for the scipy algorithm the precision must be 53"
    890         import scipy.special
    891         ans = str(scipy.special.kv(float(nu),float(z)))
    892         ans = ans.replace("(","")
    893         ans = ans.replace(")","")
    894         ans = ans.replace("j","*I")
    895         ans = sage_eval(ans)
    896         return real(ans) if z in RR else ans
    897     elif algorithm == 'pari':
    898         from sage.libs.pari.all import pari
    899         try:
    900             R = RealField(prec)
    901             nu = R(nu)
    902             z = R(z)
    903         except TypeError:
    904             C = ComplexField(prec)
    905             nu = C(nu)
    906             z = C(z)
    907             K = C
    908         K = z.parent()
    909         return K(pari(nu).besselk(z, precision=prec))
    910     elif algorithm == 'maxima':
    911         raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
    912     else:
    913         raise ValueError, "unknown algorithm '%s'"%algorithm
    914    
    915 
    916 def bessel_Y(nu,z,algorithm="maxima", prec=53):
    917     r"""
    918     Implements the "Y-Bessel function", or "Bessel function of the 2nd
    919     kind", with index (or "order") nu and argument z.
    920    
    921     .. note::
    922 
    923        Currently only prec=53 is supported.
    924    
    925     Defn::
    926    
    927                     cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
    928                    -------------------------------------------------
    929                                      sin(nu*pi)
    930    
    931     if nu is not an integer and by taking a limit otherwise.
    932    
    933     Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
    934    
    935     This is computed using Maxima by default.
    936    
    937     EXAMPLES::
    938    
    939         sage: bessel_Y(2,1.1,"scipy")
    940         -1.4314714939...
    941         sage: bessel_Y(2,1.1)   
    942         -1.4314714939590...
    943         sage: bessel_Y(3.001,2.1)
    944         -1.0299574976424...
    945 
    946     TESTS::
    947 
    948         sage: bessel_Y(2,1.1, algorithm="pari")
    949         Traceback (most recent call last):
    950         ...
    951         NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
    952     """
    953     if algorithm=="scipy":
    954         if prec != 53:
    955             raise ValueError, "for the scipy algorithm the precision must be 53"
    956         import scipy.special
    957         ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
    958         ans = ans.replace("(","")
    959         ans = ans.replace(")","")
    960         ans = ans.replace("j","*I")
    961         ans = sage_eval(ans)
    962         return real(ans) if z in RR else ans
    963     elif algorithm == "maxima":
    964         if prec != 53:
    965             raise ValueError, "for the maxima algorithm the precision must be 53"
    966         return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    967     elif algorithm == "pari":
    968         raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
    969     else:
    970         raise ValueError, "unknown algorithm '%s'"%algorithm
    971    
    972 class Bessel():
    973     """
    974     A class implementing the I, J, K, and Y Bessel functions.
    975    
    976     EXAMPLES::
    977    
    978         sage: g = Bessel(2); g
    979         J_{2}
    980         sage: print g
    981         J-Bessel function of order 2
    982         sage: g.plot(0,10)
    983 
    984     ::
    985 
    986         sage: Bessel(2, typ='I')(pi)
    987         2.61849485263445
    988         sage: Bessel(2, typ='J')(pi)
    989         0.485433932631509
    990         sage: Bessel(2, typ='K')(pi)
    991         0.0510986902537926
    992         sage: Bessel(2, typ='Y')(pi)
    993         -0.0999007139289404
    994     """
    995     def __init__(self, nu, typ = "J", algorithm = None, prec = 53):
    996         """
    997         Initializes new instance of the Bessel class.
    998 
    999         INPUT:
    1000 
    1001          - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'
    1002            or 'Y'.
    1003 
    1004          - ``algorithm`` -- (default: maxima for type Y, pari for other types)
    1005            algorithm to use to compute the Bessel function: 'pari', 'maxima' or
    1006            'scipy'.  Note that type K is not implemented in Maxima and type Y
    1007            is not implemented in PARI.
    1008 
    1009          - ``prec`` -- (default: 53) precision in bits of the Bessel function.
    1010            Only supported for the PARI algorithm.
    1011 
    1012         EXAMPLES::
    1013 
    1014             sage: g = Bessel(2); g
    1015             J_{2}
    1016             sage: Bessel(1,'I')
    1017             I_{1}
    1018             sage: Bessel(6, prec=120)(pi)
    1019             0.014545966982505560573660369604001804
    1020             sage: Bessel(6, algorithm="pari")(pi)
    1021             0.0145459669825056
    1022 
    1023         For the Bessel J-function, Maxima returns a symbolic result.  For
    1024         types I and Y, we always get a numeric result::
    1025 
    1026             sage: b = Bessel(6, algorithm="maxima")(pi); b
    1027             bessel_j(6, pi)
    1028             sage: b.n(53)
    1029             0.0145459669825056
    1030             sage: Bessel(6, typ='I', algorithm="maxima")(pi)
    1031             0.0294619840059568
    1032             sage: Bessel(6, typ='Y', algorithm="maxima")(pi)
    1033             -4.33932818939038
    1034 
    1035         SciPy usually gives less precise results::
    1036 
    1037             sage: Bessel(6, algorithm="scipy")(pi)
    1038             0.0145459669825000...
    1039 
    1040         TESTS::
    1041 
    1042             sage: Bessel(1,'Z')
    1043             Traceback (most recent call last):
    1044             ...
    1045             ValueError: typ must be one of I, J, K, Y
    1046         """
    1047         if not (typ in ['I', 'J', 'K', 'Y']):
    1048             raise ValueError, "typ must be one of I, J, K, Y"
    1049 
    1050         # Did the user ask for the default algorithm?
    1051         if algorithm is None:
    1052             if typ == 'Y':
    1053                 algorithm = 'maxima'
    1054             else:
    1055                 algorithm = 'pari'
    1056 
    1057         self._system = algorithm
    1058         self._order = nu
    1059         self._type = typ
    1060         prec = int(prec)
    1061         if prec < 0:
    1062             raise ValueError, "prec must be a positive integer"
    1063         self._prec = int(prec)
    1064 
    1065     def __str__(self):
    1066         """
    1067         Returns a string representation of this Bessel object.
    1068 
    1069         TEST::
    1070 
    1071             sage: a = Bessel(1,'I')
    1072             sage: str(a)
    1073             'I-Bessel function of order 1'
    1074         """
    1075         return self.type()+"-Bessel function of order "+str(self.order())
    1076    
    1077     def __repr__(self):
    1078         """
    1079         Returns a string representation of this Bessel object.
    1080 
    1081         TESTS::
    1082 
    1083             sage: Bessel(1,'I')
    1084             I_{1}
    1085         """
    1086         return self.type()+"_{"+str(self.order())+"}"
    1087    
    1088     def type(self):
    1089         """
    1090         Returns the type of this Bessel object.
    1091 
    1092         TEST::
    1093 
    1094             sage: a = Bessel(3,'K')
    1095             sage: a.type()
    1096             'K'
    1097         """
    1098         return self._type
    1099    
    1100     def prec(self):
    1101         """
    1102         Returns the precision (in number of bits) used to represent this
    1103         Bessel function.
    1104 
    1105         TESTS::
    1106 
    1107             sage: a = Bessel(3,'K')
    1108             sage: a.prec()
    1109             53
    1110             sage: B = Bessel(20,prec=100); B
    1111             J_{20}
    1112             sage: B.prec()
    1113             100
    1114         """
    1115         return self._prec
    1116 
    1117     def order(self):
    1118         """
    1119         Returns the order of this Bessel function.
    1120 
    1121         TEST::
    1122 
    1123             sage: a = Bessel(3,'K')
    1124             sage: a.order()
    1125             3
    1126         """
    1127         return self._order
    1128 
    1129     def system(self):
    1130         """
    1131         Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with
    1132         this Bessel function.
    1133 
    1134         TESTS::
    1135 
    1136             sage: Bessel(20,algorithm='maxima').system()
    1137             'maxima'
    1138             sage: Bessel(20,prec=100).system()
    1139             'pari'
    1140         """
    1141         return self._system
    1142 
    1143     def __call__(self,z):
    1144         """
    1145         Implements evaluation of all the Bessel functions directly
    1146         from the Bessel class. This essentially allows a Bessel object to
    1147         behave like a function that can be invoked.
    1148 
    1149         TESTS::
    1150 
    1151             sage: Bessel(3,'K')(5.0)
    1152             0.00829176841523093
    1153             sage: Bessel(20,algorithm='maxima')(5.0)
    1154             2.77033005213e-11
    1155             sage: Bessel(20,prec=100)(5.0101010101010101)
    1156             2.8809188227195382093062257967e-11
    1157             sage: B = Bessel(2,'Y',algorithm='scipy',prec=50)
    1158             sage: B(2.0)
    1159             Traceback (most recent call last):
    1160             ...
    1161             ValueError: for the scipy algorithm the precision must be 53
    1162         """
    1163         nu = self.order()
    1164         t = self.type()
    1165         s = self.system()
    1166         p = self.prec()
    1167         if t == "I":
    1168             return bessel_I(nu,z,algorithm=s,prec=p)
    1169         if t == "J":
    1170             return bessel_J(nu,z,algorithm=s,prec=p)
    1171         if t == "K":
    1172             return bessel_K(nu,z,algorithm=s,prec=p)
    1173         if t == "Y":
    1174             return bessel_Y(nu,z,algorithm=s,prec=p)
    1175        
    1176     def plot(self,a,b):
    1177         """
    1178         Enables easy plotting of all the Bessel functions directly
    1179         from the Bessel class.
    1180 
    1181         TESTS::
    1182 
    1183             sage: plot(Bessel(2),3,4)
    1184             sage: Bessel(2).plot(3,4)
    1185             sage: P = Bessel(2,'I').plot(1,5)
    1186             sage: P += Bessel(2,'J').plot(1,5)
    1187             sage: P += Bessel(2,'K').plot(1,5)
    1188             sage: P += Bessel(2,'Y').plot(1,5)
    1189             sage: show(P)
    1190         """
    1191         nu = self.order()
    1192         s = self.system()
    1193         t = self.type()
    1194         if t == "I":
    1195             f = lambda z: bessel_I(nu,z,s) 
    1196             P = plot(f,a,b)
    1197         if t == "J":
    1198             f = lambda z: bessel_J(nu,z,s)
    1199             P = plot(f,a,b)
    1200         if t == "K":
    1201             f = lambda z: bessel_K(nu,z,s)
    1202             P = plot(f,a,b)
    1203         if t == "Y":
    1204             f = lambda z: bessel_Y(nu,z,s)
    1205             P = plot(f,a,b)
    1206         return P
    1207    
    1208606def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53):
    1209607    r"""
    1210608    Default is a wrap of PARI's hyperu(alpha,beta,x) function.