Ticket #4102: trac_symbolic_bessel_v6.patch

File trac_symbolic_bessel_v6.patch, 50.8 KB (added by benjaminfjones, 8 years ago)

Copy of trac_symbolic_bessel_v5, minus one doctest block

  • doc/en/reference/functions/index.rst

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1363120860 25200
    # Node ID b0469e15629643dfc405ec3f2f60705b9c1c1742
    # Parent  1077314f416653b28e199c382667a1f11e444bdd
    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 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.rings.all import RR, Integer
     176from sage.misc.latex import latex
     177from sage.symbolic.function import BuiltinFunction, is_inexact
     178from sage.symbolic.expression import Expression
     179from sage.structure.coerce import parent
     180import sage.structure.element
     181from sage.libs.mpmath import utils as mpmath_utils
     182from sage.functions.other import sqrt
     183from sage.functions.log import exp
     184from sage.functions.hyperbolic import sinh, cosh
     185from sage.symbolic.constants import pi
     186
     187
     188class Function_Bessel_J(BuiltinFunction):
     189    r"""
     190    The Bessel J Function, denoted by bessel_J(`\nu`, x) or `J_\nu(x)`.
     191    As a Taylor series about `x=0` it is equal to:
     192
     193    .. math::
     194
     195        J_\nu(x) = \sum_{k=0}^\infty \frac{(-1)^k}{k! \Gamma(k+\nu+1)}
     196        (\frac{x}{2})^{2k+\nu}
     197
     198    The parameter `\nu` is called the order and may be any real or
     199    complex number, however the integer and half-integer values are most
     200    common. It is defined for all complex numbers `x` when `\nu`
     201    is an integer or greater than zero and it diverges as `x \to 0`
     202    for negative non-integer values of `\nu`.
     203
     204    For integer orders `\nu = n` there is an integral representation:
     205
     206    .. math::
     207
     208        J_n(x) = \frac{1}{\pi} \int_0^\pi \cos(n t - x \sin(t)) \; dt
     209
     210    This function also arises as a special case of the hypergeometric
     211    function `{}_0F_1`:
     212
     213    .. math::
     214
     215        J_\nu(x) = \frac{x^n}{2^\nu \Gamma(\nu + 1)} {}_0F_1(\nu +
     216        1, -\frac{x^2}{4}).
     217
     218    EXAMPLES::
     219
     220        sage: bessel_J(1.0, 1.0)
     221        0.440050585744933
     222        sage: bessel_J(2, I).n(digits=30)
     223        -0.135747669767038281182852569995
     224
     225        sage: bessel_J(1, x)
     226        bessel_J(1, x)
     227        sage: n = var('n')
     228        sage: bessel_J(n, x)
     229        bessel_J(n, x)
     230
     231    Examples of symbolic manipulation::
     232
     233        sage: a = bessel_J(pi, bessel_J(1, I)); a
     234        bessel_J(pi, bessel_J(1, I))
     235        sage: N(a, digits=20)
     236        0.00059023706363796717363 - 0.0026098820470081958110*I
     237
     238        sage: f = bessel_J(2, x)
     239        sage: f.diff(x)
     240        1/2*bessel_J(1, x) - 1/2*bessel_J(3, x)
     241
     242    Comparison to a well-known integral representation of `J_1(1)`::
     243
     244        sage: A = numerical_integral(1/pi*cos(x - sin(x)), 0, pi)
     245        sage: A[0]
     246        0.44005058574493355
     247        sage: bessel_J(1.0, 1.0) - A[0] < 1e-15
     248        True
     249
     250    Currently, integration is not supported (directly) since we cannot
     251    yet convert hypergeometric functions to and from Maxima::
     252
     253        sage: f = bessel_J(2, x)
     254        sage: f.integrate(x)
     255        Traceback (most recent call last):
     256        ...
     257        TypeError: cannot coerce arguments: no canonical coercion from <type 'list'> to Symbolic Ring
     258
     259        sage: m = maxima(bessel_J(2, x))
     260        sage: m.integrate(x)
     261        hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24
     262
     263    Visualization::
     264
     265        sage: plot(bessel_J(1,x), (x,0,5), color='blue')
     266        sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5)) # long time
     267
     268    ALGORITHM:
     269
     270        Numerical evaluation is handled by the mpmath library. Symbolics are
     271        handled by a combination of Maxima and Sage (Ginac/Pynac).
     272    """
     273    def __init__(self):
     274        """
     275        See the docstring for :meth:`Function_Bessel_J`.
     276
     277        EXAMPLES::
     278
     279            sage: sage.functions.bessel.Function_Bessel_J()
     280            bessel_J
     281        """
     282        BuiltinFunction.__init__(self, "bessel_J", nargs=2,
     283                                 conversions=dict(maxima='bessel_j',
     284                                                  mathematica='BesselJ'))
     285
     286    def _eval_(self, n, x):
     287        """
     288        EXAMPLES::
     289
     290            sage: a, b = var('a, b')
     291            sage: bessel_J(a, b)
     292            bessel_J(a, b)
     293            sage: bessel_J(1.0, 1.0)
     294            0.440050585744933
     295        """
     296        if (not isinstance(n, Expression) and
     297                not isinstance(x, Expression) and
     298                (is_inexact(n) or is_inexact(x))):
     299            coercion_model = sage.structure.element.get_coercion_model()
     300            n, x = coercion_model.canonical_coercion(n, x)
     301            return self._evalf_(n, x, parent(n))
     302
     303        return None
     304
     305    def _evalf_(self, n, x, parent=None):
     306        """
     307        EXAMPLES::
     308
     309            sage: bessel_J(0.0, 1.0)
     310            0.765197686557966
     311            sage: bessel_J(0, 1).n(digits=20)
     312            0.76519768655796655145
     313        """
     314        import mpmath
     315        return mpmath_utils.call(mpmath.besselj, n, x, parent=parent)
     316
     317    def _derivative_(self, n, x, diff_param=None):
     318        """
     319        Return the derivative of the Bessel J function.
     320
     321        EXAMPLES::
     322
     323            sage: f(z) = bessel_J(10, z)
     324            sage: derivative(f, z)
     325            z |--> 1/2*bessel_J(9, z) - 1/2*bessel_J(11, z)
     326        """
     327        return (bessel_J(n-1, x) - bessel_J(n+1, x))/Integer(2)
     328
     329    def _print_latex_(self, n, z):
     330        """
     331        Custom _print_latex_ method.
     332
     333        EXAMPLES::
     334
     335            sage: latex(bessel_J(1, x))
     336            \operatorname{J_{1}}(x)
     337        """
     338        return r"\operatorname{J_{%s}}(%s)" % (latex(n), latex(z))
     339
     340bessel_J = Function_Bessel_J()
     341
     342
     343class Function_Bessel_Y(BuiltinFunction):
     344    r"""
     345    The Bessel Y functions, also known as the Bessel functions of the second
     346    kind, Weber functions, or Neumann functions.
     347
     348    `Y_\nu(z)` is a holomorphic function of `z` on the complex plane,
     349    cut along the negative real axis. It is singular at `z = 0`. When `z`
     350    is fixed, `Y_\nu(z)` is an entire function of the order `\nu`.
     351
     352    DEFINITION:
     353
     354    .. math::
     355
     356        Y_n(z) = \frac{J_\nu(z) \cos(\nu z) -
     357        J_{-\nu}(z)}{\sin(\nu z)}
     358
     359    Its derivative with respect to `z` is:
     360
     361    .. math::
     362
     363        \frac{d}{dz} Y_n(z) = \frac{1}{z^n} \left(z^n Y_{n-1}(z) - n z^{n-1}
     364        Y_n(z) \right)
     365
     366    EXAMPLES::
     367
     368        sage: bessel_Y(1, x)
     369        bessel_Y(1, x)
     370        sage: bessel_Y(1.0, 1.0)
     371        -0.781212821300289
     372        sage: n = var('n')
     373        sage: bessel_Y(n, x)
     374        bessel_Y(n, x)
     375        sage: bessel_Y(2, I).n()
     376        1.03440456978312 - 0.135747669767038*I
     377        sage: bessel_Y(0, 0).n()
     378        -infinity
     379
     380    Examples of symbolic manipulation::
     381
     382        sage: a = bessel_Y(pi, bessel_Y(1, I)); a
     383        bessel_Y(pi, bessel_Y(1, I))
     384        sage: N(a, digits=20)
     385        4.2059146571791095708 + 21.307914215321993526*I
     386
     387        sage: f = bessel_Y(2, x)
     388        sage: f.diff(x)
     389        1/2*bessel_Y(1, x) - 1/2*bessel_Y(3, x)
     390
     391    High precision and complex valued inputs (see :trac:`4230`)::
     392
     393        sage: bessel_Y(0, 1).n(128)
     394        0.088256964215676957982926766023515162828
     395        sage: bessel_Y(0, RealField(200)(1))
     396        0.088256964215676957982926766023515162827817523090675546711044
     397        sage: bessel_Y(0, ComplexField(200)(0.5+I))
     398        0.077763160184438051408593468823822434235010300228009867784073 + 1.0142336049916069152644677682828326441579314239591288411739*I
     399
     400    Visualization::
     401
     402        sage: plot(bessel_Y(1,x), (x,0,5), color='blue')
     403        sage: complex_plot(bessel_Y(1, x), (-5, 5), (-5, 5)) # long time
     404
     405    ALGORITHM:
     406
     407        Numerical evaluation is handled by the mpmath library. Symbolics are
     408        handled by a combination of Maxima and Sage (Ginac/Pynac).
     409    """
     410    def __init__(self):
     411        """
     412        See the docstring for :meth:`Function_Bessel_Y`.
     413
     414        EXAMPLES::
     415
     416            sage: sage.functions.bessel.Function_Bessel_Y()(0, x)
     417            bessel_Y(0, x)
     418        """
     419        BuiltinFunction.__init__(self, "bessel_Y", nargs=2,
     420                                 conversions=dict(maxima='bessel_y',
     421                                                  mathematica='BesselY'))
     422
     423    def _eval_(self, n, x):
     424        """
     425        EXAMPLES::
     426
     427            sage: a,b = var('a, b')
     428            sage: bessel_Y(a, b)
     429            bessel_Y(a, b)
     430            sage: bessel_Y(0, 1).n(128)
     431            0.088256964215676957982926766023515162828
     432        """
     433        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     434                (is_inexact(n) or is_inexact(x))):
     435            coercion_model = sage.structure.element.get_coercion_model()
     436            n, x = coercion_model.canonical_coercion(n, x)
     437            return self._evalf_(n, x, parent(n))
     438
     439        return None  # leaves the expression unevaluated
     440
     441    def _evalf_(self, n, x, parent=None):
     442        """
     443        EXAMPLES::
     444
     445            sage: bessel_Y(1.0+2*I, 3.0+4*I)
     446            0.699410324467538 + 0.228917940896421*I
     447            sage: bessel_Y(0, 1).n(256)
     448            0.08825696421567695798292676602351516282781752309067554671104384761199978932351
     449        """
     450        import mpmath
     451        return mpmath_utils.call(mpmath.bessely, n, x, parent=parent)
     452
     453    def _derivative_(self, n, x, diff_param=None):
     454        """
     455        Return the derivative of the Bessel Y function.
     456
     457        EXAMPLES::
     458
     459            sage: f(x) = bessel_Y(10, x)
     460            sage: derivative(f, x)
     461            x |--> 1/2*bessel_Y(9, x) - 1/2*bessel_Y(11, x)
     462        """
     463        return (bessel_Y(n-1, x) - bessel_Y(n+1, x))/Integer(2)
     464
     465    def _print_latex_(self, n, z):
     466        """
     467        Custom _print_latex_ method.
     468
     469        EXAMPLES::
     470
     471            sage: latex(bessel_Y(1, x))
     472            \operatorname{Y_{1}}(x)
     473        """
     474        return r"\operatorname{Y_{%s}}(%s)" % (latex(n), latex(z))
     475
     476bessel_Y = Function_Bessel_Y()
     477
     478
     479class Function_Bessel_I(BuiltinFunction):
     480    r"""
     481    The Bessel I function, or the Modified Bessel Function of the First Kind.
     482
     483    DEFINITION:
     484
     485    .. math::
     486
     487        I_\nu(x) = i^{-\nu} J_\nu(ix)
     488
     489    EXAMPLES::
     490
     491        sage: bessel_I(1, x)
     492        bessel_I(1, x)
     493        sage: bessel_I(1.0, 1.0)
     494        0.565159103992485
     495        sage: n = var('n')
     496        sage: bessel_I(n, x)
     497        bessel_I(n, x)
     498        sage: bessel_I(2, I).n()
     499        -0.114903484931900
     500
     501    Examples of symbolic manipulation::
     502
     503        sage: a = bessel_I(pi, bessel_I(1, I))
     504        sage: N(a, digits=20)
     505        0.00026073272117205890528 - 0.0011528954889080572266*I
     506
     507        sage: f = bessel_I(2, x)
     508        sage: f.diff(x)
     509        1/2*bessel_I(1, x) + 1/2*bessel_I(3, x)
     510
     511    Special identities that bessel_I satisfies::
     512
     513        sage: bessel_I(1/2, x)
     514        sqrt(1/(pi*x))*sqrt(2)*sinh(x)
     515        sage: eq = bessel_I(1/2, x) == bessel_I(0.5, x)
     516        sage: eq.test_relation()
     517        True
     518        sage: bessel_I(-1/2, x)
     519        sqrt(1/(pi*x))*sqrt(2)*cosh(x)
     520        sage: eq = bessel_I(-1/2, x) == bessel_I(-0.5, x)
     521        sage: eq.test_relation()
     522        True
     523
     524    Examples of asymptotic behavior::
     525
     526        sage: limit(bessel_I(0, x), x=oo)
     527        +Infinity
     528        sage: limit(bessel_I(0, x), x=0)
     529        1
     530
     531    High precision and complex valued inputs::
     532
     533        sage: bessel_I(0, 1).n(128)
     534        1.2660658777520083355982446252147175376
     535        sage: bessel_I(0, RealField(200)(1))
     536        1.2660658777520083355982446252147175376076703113549622068081
     537        sage: bessel_I(0, ComplexField(200)(0.5+I))
     538        0.80644357583493619472428518415019222845373366024179916785502 + 0.22686958987911161141397453401487525043310874687430711021434*I
     539
     540    Visualization::
     541
     542        sage: plot(bessel_I(1,x), (x,0,5), color='blue')
     543        sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5)) # long time
     544
     545    ALGORITHM:
     546
     547        Numerical evaluation is handled by the mpmath library. Symbolics are
     548        handled by a combination of Maxima and Sage (Ginac/Pynac).
     549    """
     550    def __init__(self):
     551        """
     552        See the docstring for :meth:`Function_Bessel_I`.
     553
     554        EXAMPLES::
     555
     556            sage: bessel_I(1,x)
     557            bessel_I(1, x)
     558        """
     559        BuiltinFunction.__init__(self, "bessel_I", nargs=2,
     560                                 conversions=dict(maxima='bessel_i',
     561                                                  mathematica='BesselI'))
     562
     563    def _eval_(self, n, x):
     564        """
     565        EXAMPLES::
     566
     567            sage: y=var('y')
     568            sage: bessel_I(y,x)
     569            bessel_I(y, x)
     570            sage: bessel_I(0.0, 1.0)
     571            1.26606587775201
     572            sage: bessel_I(1/2, 1)
     573            sqrt(2)*sinh(1)/sqrt(pi)
     574            sage: bessel_I(-1/2, pi)
     575            sqrt(2)*cosh(pi)/pi
     576        """
     577        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     578                (is_inexact(n) or is_inexact(x))):
     579            coercion_model = sage.structure.element.get_coercion_model()
     580            n, x = coercion_model.canonical_coercion(n, x)
     581            return self._evalf_(n, x, parent(n))
     582
     583        # special identities
     584        if n == Integer(1)/Integer(2):
     585            return sqrt(2/(pi*x))*sinh(x)
     586        elif n == -Integer(1)/Integer(2):
     587            return sqrt(2/(pi*x))*cosh(x)
     588
     589        return None  # leaves the expression unevaluated
     590
     591    def _evalf_(self, n, x, parent=None):
     592        """
     593        EXAMPLES::
     594
     595            sage: bessel_I(1,3).n(digits=20)
     596            3.9533702174026093965
     597        """
     598        import mpmath
     599        return mpmath_utils.call(mpmath.besseli, n, x, parent=parent)
     600
     601    def _derivative_(self, n, x, diff_param=None):
     602        """
     603        Return the derivative of the Bessel I function `I_n(x)` with repect
     604        to `x`.
     605
     606        EXAMPLES::
     607
     608            sage: f(z) = bessel_I(10, x)
     609            sage: derivative(f, x)
     610            z |--> 1/2*bessel_I(9, x) + 1/2*bessel_I(11, x)
     611        """
     612        return (bessel_I(n-1, x) + bessel_I(n+1, x))/Integer(2)
     613
     614    def _print_latex_(self, n, z):
     615        """
     616        Custom _print_latex_ method.
     617
     618        EXAMPLES::
     619
     620            sage: latex(bessel_I(1, x))
     621            \operatorname{I_{1}}(x)
     622        """
     623        return r"\operatorname{I_{%s}}(%s)" % (latex(n), latex(z))
     624
     625bessel_I = Function_Bessel_I()
     626
     627
     628class Function_Bessel_K(BuiltinFunction):
     629    r"""
     630    The Bessel K function, or the modified Bessel function of the second kind.
     631
     632    DEFINITION:
     633
     634    .. math::
     635
     636        K_\nu(x) = \frac{\pi}{2} \frac{I_{-\nu}(x)-I_\nu(x)}{\sin(\nu \pi)}
     637
     638    EXAMPLES::
     639
     640        sage: bessel_K(1, x)
     641        bessel_K(1, x)
     642        sage: bessel_K(1.0, 1.0)
     643        0.601907230197235
     644        sage: n = var('n')
     645        sage: bessel_K(n, x)
     646        bessel_K(n, x)
     647        sage: bessel_K(2, I).n()
     648        -2.59288617549120 + 0.180489972066962*I
     649
     650    Examples of symbolic manipulation::
     651
     652        sage: a = bessel_K(pi, bessel_K(1, I)); a
     653        bessel_K(pi, bessel_K(1, I))
     654        sage: N(a, digits=20)
     655        3.8507583115005220157 + 0.068528298579883425792*I
     656
     657        sage: f = bessel_K(2, x)
     658        sage: f.diff(x)
     659        1/2*bessel_K(1, x) + 1/2*bessel_K(3, x)
     660
     661        sage: bessel_K(1/2, x)
     662        bessel_K(1/2, x)
     663        sage: bessel_K(1/2, -1)
     664        bessel_K(1/2, -1)
     665        sage: bessel_K(1/2, 1)
     666        sqrt(pi)*sqrt(1/2)*e^(-1)
     667
     668    Examples of asymptotic behavior::
     669
     670        sage: bessel_K(0, 0.0)
     671        +infinity
     672        sage: limit(bessel_K(0, x), x=0)
     673        +Infinity
     674        sage: limit(bessel_K(0, x), x=oo)
     675        0
     676
     677    High precision and complex valued inputs::
     678
     679        sage: bessel_K(0, 1).n(128)
     680        0.42102443824070833333562737921260903614
     681        sage: bessel_K(0, RealField(200)(1))
     682        0.42102443824070833333562737921260903613621974822666047229897
     683        sage: bessel_K(0, ComplexField(200)(0.5+I))
     684        0.058365979093103864080375311643360048144715516692187818271179 - 0.67645499731334483535184142196073004335768129348518210260256*I
     685
     686    Visualization::
     687
     688        sage: plot(bessel_K(1,x), (x,0,5), color='blue')
     689        sage: complex_plot(bessel_K(1, x), (-5, 5), (-5, 5)) # long time
     690
     691    ALGORITHM:
     692
     693        Numerical evaluation is handled by the mpmath library. Symbolics are
     694        handled by a combination of Maxima and Sage (Ginac/Pynac).
     695
     696    TESTS:
     697
     698    Verify that :trac:`3426` is fixed:
     699
     700    The Bessel K function can be evaluated numerically at complex orders::
     701
     702        sage: bessel_K(10 * I, 10).n()
     703        9.82415743819925e-8
     704
     705    For a fixed imaginary order and increasing, real, second component the
     706    value of Bessel K is exponentially decaying::
     707
     708        sage: for x in [10, 20, 50, 100, 200]: print bessel_K(5*I, x).n()
     709        5.27812176514912e-6
     710        3.11005908421801e-10
     711        2.66182488515423e-23 - 8.59622057747552e-58*I
     712        4.11189776828337e-45 - 1.01494840019482e-80*I
     713        1.15159692553603e-88 - 6.75787862113718e-125*I
     714    """
     715    def __init__(self):
     716        """
     717        See the docstring for :meth:`Function_Bessel_K`.
     718
     719        EXAMPLES::
     720
     721            sage: sage.functions.bessel.Function_Bessel_K()
     722            bessel_K
     723        """
     724        BuiltinFunction.__init__(self, "bessel_K", nargs=2,
     725                                 conversions=dict(maxima='bessel_k',
     726                                                  mathematica='BesselK'))
     727
     728    def _eval_(self, n, x):
     729        """
     730        EXAMPLES::
     731
     732            sage: bessel_K(1,0)
     733            bessel_K(1, 0)
     734            sage: bessel_K(1.0, 0.0)
     735            +infinity
     736            sage: bessel_K(-1, 1).n(128)
     737            0.60190723019723457473754000153561733926
     738        """
     739        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
     740                (is_inexact(n) or is_inexact(x))):
     741            coercion_model = sage.structure.element.get_coercion_model()
     742            n, x = coercion_model.canonical_coercion(n, x)
     743            return self._evalf_(n, x, parent(n))
     744
     745        # special identity
     746        if n == Integer(1)/Integer(2) and x > 0:
     747            return sqrt(pi/2)*exp(-x)*x**(Integer(1)/Integer(2))
     748
     749        return None  # leaves the expression unevaluated
     750
     751    def _evalf_(self, n, x, parent=None):
     752        """
     753        EXAMPLES::
     754
     755            sage: bessel_K(0.0, 1.0)
     756            0.421024438240708
     757            sage: bessel_K(0, RealField(128)(1))
     758            0.42102443824070833333562737921260903614
     759        """
     760        import mpmath
     761        return mpmath_utils.call(mpmath.besselk, n, x, parent=parent)
     762
     763    def _derivative_(self, n, x, diff_param=None):
     764        """
     765        Return the derivative of the Bessel K function.
     766
     767        EXAMPLES::
     768
     769            sage: f(x) = bessel_K(10, x)
     770            sage: derivative(f, x)
     771            x |--> 1/2*bessel_K(9, x) + 1/2*bessel_K(11, x)
     772        """
     773        return (bessel_K(n-1, x) + bessel_K(n+1, x))/Integer(2)
     774
     775    def _print_latex_(self, n, z):
     776        """
     777        Custom _print_latex_ method.
     778
     779        EXAMPLES::
     780
     781            sage: latex(bessel_K(1, x))
     782            \operatorname{K_{1}}(x)
     783        """
     784        return r"\operatorname{K_{%s}}(%s)" % (latex(n), latex(z))
     785
     786bessel_K = Function_Bessel_K()
     787
     788
     789# dictionary used in Bessel
     790bessel_type_dict = {'I': bessel_I, 'J': bessel_J, 'K': bessel_K, 'Y': bessel_Y}
     791
     792
     793def Bessel(*args, **kwds):
     794    """
     795    A function factory that produces symbolic I, J, K, and Y Bessel functions.
     796    There are several ways to call this function:
     797
     798        - ``Bessel(order, type)``
     799        - ``Bessel(order)`` -- type defaults to 'J'
     800        - ``Bessel(order, typ=T)``
     801        - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter
     802          function
     803        - ``Bessel()`` -- order is unspecified, type is 'J'
     804
     805    where `order` can be any integer and T must be one of the strings 'I', 'J',
     806    'K', or 'Y'.
     807
     808    See the EXAMPLES below.
     809
     810    EXAMPLES:
     811
     812    Construction of Bessel functions with various orders and types::
     813
     814        sage: Bessel()
     815        bessel_J
     816        sage: Bessel(1)(x)
     817        bessel_J(1, x)
     818        sage: Bessel(1, 'Y')(x)
     819        bessel_Y(1, x)
     820        sage: Bessel(-2, 'Y')(x)
     821        bessel_Y(-2, x)
     822        sage: Bessel(typ='K')
     823        bessel_K
     824        sage: Bessel(0, typ='I')(x)
     825        bessel_I(0, x)
     826
     827    Evaluation::
     828
     829        sage: f = Bessel(1)
     830        sage: f(3.0)
     831        0.339058958525936
     832        sage: f(3)
     833        bessel_J(1, 3)
     834        sage: f(3).n(digits=50)
     835        0.33905895852593645892551459720647889697308041819801
     836
     837        sage: g = Bessel(typ='J')
     838        sage: g(1,3)
     839        bessel_J(1, 3)
     840        sage: g(2, 3+I).n()
     841        0.634160370148554 + 0.0253384000032695*I
     842        sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15
     843        True
     844
     845    Symbolic calculus::
     846
     847        sage: f(x) = Bessel(0, 'J')(x)
     848        sage: derivative(f, x)
     849        x |--> 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x)
     850        sage: derivative(f, x, x)
     851        x |--> 1/4*bessel_J(-2, x) - 1/2*bessel_J(0, x) + 1/4*bessel_J(2, x)
     852
     853    Verify that `J_0` satisfies Bessel's differential equation numerically
     854    using the ``test_relation()`` method::
     855
     856        sage: y = bessel_J(0, x)
     857        sage: diffeq = x^2*derivative(y,x,x) + x*derivative(y,x) + x^2*y == 0
     858        sage: diffeq.test_relation(proof=False)
     859        True
     860
     861    Conversion to other systems::
     862
     863        sage: x,y = var('x,y')
     864        sage: f = maxima(Bessel(typ='K')(x,y))
     865        sage: f.derivative('x')
     866        %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)
     867        sage: f.derivative('y')
     868        -(bessel_k(x+1,y)+bessel_k(x-1,y))/2
     869
     870    Plotting::
     871
     872        sage: f(x) = Bessel(0)(x); f
     873        x |--> bessel_J(0, x)
     874        sage: plot(f, (x, 1, 10))
     875
     876        sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10)
     877
     878        sage: G = Graphics()
     879        sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ])
     880        sage: show(G)
     881
     882    A recreation of Abramowitz and Stegun Figure 9.1::
     883
     884        sage: G  = plot(Bessel(0, 'J'), 0, 15, color='black')
     885        sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black')
     886        sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted')
     887        sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted')
     888        sage: show(G, ymin=-1, ymax=1)
     889
     890    """
     891    # Determine the order and type of function from the arguments and keywords.
     892    # These are recored in local variables: _type, _order, _system, _nargs.
     893    _type = None
     894    if len(args) == 0:    # no order specified
     895        _order = None
     896        _nargs = 2
     897    elif len(args) == 1:  # order is specified
     898        _order = args[0]
     899        _nargs = 1
     900    elif len(args) == 2:  # both order and type are positional arguments
     901        _order = args[0]
     902        _type = args[1]
     903        _nargs = 1
     904    else:
     905        raise TypeError("at most two positional arguments may be specified, "
     906                        + "see the docstring for Bessel")
     907    # check for type inconsistency
     908    if _type is not None and 'typ' in kwds and _type != kwds['typ']:
     909        raise ValueError("inconsistent types given")
     910    # record the function type
     911    if _type is None:
     912        if 'typ' in kwds:
     913            _type = kwds['typ']
     914        else:
     915            _type = 'J'
     916    if not (_type in ['I', 'J', 'K', 'Y']):
     917        raise ValueError("type must be one of I, J, K, Y")
     918    # record the numerical evaluation system
     919    if 'algorithm' in kwds:
     920        _system = kwds['algorithm']
     921    else:
     922        _system = 'mpmath'
     923
     924    # return the function
     925    # TODO remove assertions
     926    assert _type in ['I', 'J', 'K', 'Y']
     927    assert _order is None or _order in RR
     928    assert _nargs == 1 or _nargs == 2
     929    assert _system == 'mpmath'
     930
     931    # TODO what to do with _system?
     932    _f = bessel_type_dict[_type]
     933    if _nargs == 1:
     934        return lambda x: _f(_order, x)
     935    else:
     936        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
    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.