Ticket #4102: trac_symbolic_bessel_v3.patch

File trac_symbolic_bessel_v3.patch, 56.5 KB (added by benjaminfjones, 9 years ago)

more work in progress, v3

  • doc/en/reference/functions.rst 

    # HG changeset patch
# User Benjamin Jones <benjaminfjones@gmail.com>
# Date 1356654134 28800
# Node ID eb9efa4564e2b4f095c2abc0796e81fed6ec4510
# Parent  3df9be5aadd2e4159e69bdf911d38784cb6fa130
Trac 4102: Implement symbolic Bessel functions

diff --git a/doc/en/reference/functions.rst b/doc/en/reference/functions.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.
     6
     7-  Bessel functions, first defined by the Swiss mathematician
     8   Daniel Bernoulli and named after Friedrich Bessel, are canonical
     9   solutions y(x) of Bessel's differential equation:
     10
     11   .. math::
     12
     13         x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0,
     14
     15   for an arbitrary real number `\alpha` (the order).
     16
     17-  In this module, `J_\alpha` denotes the unqiue solution of Bessel's equation
     18   which is non-singular at `x = 0`. This function is known as the Bessel
     19   Function of the First Kind. This function also arrises as a special case
     20   of the hypergeometric function `{}_0F_1`:
     21
     22   .. math::
     23
     24        J_\alpha(x) = \frac{x^n}{2^\alpha \Gamma(\alpha + 1)} {}_0F_1(\alpha + 1, -\frac{x^2}{4}).
     25
     26-  The second linearly independent solution to Bessel's equation (which is
     27   singular at `x=0`) is denoted by `Y_\alpha` and is called the Bessel Function
     28   of the Second Kind:
     29
     30   .. math::
     31
     32        Y_\alpha(x) = \frac{ J_\alpha(x) \cos(\pi \alpha) - J_{-\alpha}(x)}{\sin(\pi \alpha)}.
     33
     34-  There are also two commonly used combinations of the Bessel J and Y Functions.
     35   The Bessel I Function, or the Modified Bessel Function of the First Kind, is
     36   defined by:
     37
     38   .. math::
     39
     40       I_\alpha(x) = i^{-\alpha} J_\alpha(ix).
     41
     42   The Bessel K Function, or the Modified Bessel Function of the Second Kind,
     43   is deinfed by:
     44
     45   .. math::
     46
     47       K_\alpha(x) = \frac{\pi}{2} \cdot \frac{I_{-\alpha}(x) - I_n(x)}{\sin(\pi \alpha)}.
     48
     49   We should note here that the above formulas for Bessel Y and K functions
     50   should be understood as limits when `\alpha` is an integer.
     51
     52-  It follows from Bessel's differential equation that the derivative of
     53   `J_n(x)` with respect to `x` is:
     54
     55    ..math::
     56
     57        \frac{d}{dx} J_n(x) = \frac{1}{x^n} \left\(x^n J_{n-1}(x) - n x^{n-1} J_n(z) \right\)
     58
     59-  Another important formulation of the two linearly independent
     60   solutions to Bessel's equation are the Hankel functions
     61   `H_\alpha^{(1)}(x)` and `H_\alpha^{(2)}(x)`,
     62   defined by:
     63
     64   .. math::
     65
     66         H_\alpha^{(1)}(x) = J_\alpha(x) + i Y_\alpha(x)
     67
     68   .. math::
     69
     70         H_\alpha^{(2)}(x) = J_\alpha(x) - i Y_\alpha(x)
     71
     72   where `i` is the imaginary unit (and `J_*` and
     73   `Y_*` are the usual J- and Y-Bessel functions). These
     74   linear combinations are also known as Bessel functions of the third
     75   kind; they are two linearly independent solutions of Bessel's
     76   differential equation. They are named for Hermann Hankel.
     77
     78EXAMPLES:
     79
     80    Evaluate the Bessel J function symbolically and numerically::
     81
     82        sage: bessel_J(0, x)
     83        bessel_J(0, x)
     84        sage: bessel_J(0, 0)
     85        bessel_J(0, 0)
     86        sage: bessel_J(0, x).diff(x)
     87        1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x)
     88
     89        sage: N(bessel_J(0, 0), digits = 20)
     90        1.0000000000000000000
     91        sage: find_root(bessel_J(0,x), 0, 5)
     92        2.404825557695773
     93
     94    Plot the Bessel J function::
     95
     96        sage: f(x) = Bessel(0)(x); f
     97        x |--> bessel_J(0, x)
     98        sage: plot(f, (x, 1, 10))
     99
     100    Visualize the Bessel Y function on the complex plane::
     101
     102        sage: complex_plot(bessel_Y(0, x), (-5, 5), (-5, 5))
     103
     104    Evaluate a combination of Bessel functions::
     105
     106        sage: f(x) = bessel_J(1, x) - bessel_Y(0, x)
     107        sage: f(pi)
     108        bessel_J(1, pi) - bessel_Y(0, pi)
     109        sage: f(pi).n()
     110        -0.0437509653365599
     111        sage: f(pi).n(digits=50)
     112        -0.043750965336559909054985168023342675387737118378169
     113
     114    Symbolically solve a second order differential equation with initial
     115    conditions `y(1) = a` and `y'(1) = b` in terms of Bessel functions::
     116
     117        sage: y = function('y', x)
     118        sage: a, b = var('a, b')
     119        sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0
     120        sage: f = desolve(diffeq, y, [1, a, b]); f
     121        (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))
     122
     123    For more examples, see the docstring for :meth:`Bessel`.
     124
     125AUTHORS:
     126
     127    - Benjamin Jones (2012-12-27): initial version.
     128
     129    - Some of the documentation here has been adapted from David Joyner's
     130      original documentation of Sage's special functions module (2006).
     131
     132REFERENCES:
     133
     134    - Abramowitz and Stegun: Handbook of Mathematical Functions,
     135      http://www.math.sfu.ca/~cbm/aands/
     136
     137    - http://en.wikipedia.org/wiki/Bessel_function
     138
     139    - mpmath Library `Bessel Functions`_
     140
     141.. _`mpmath Library`: http://code.google.com/p/mpmath/
     142.. _`Bessel Functions`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html
     143
     144"""
     145
     146#*****************************************************************************
     147#       Copyright (C) 2012 Benjamin Jones <benjaminfjones@gmail.com>
     148#
     149#  Distributed under the terms of the GNU General Public License (GPL)
     150#
     151#    This code is distributed in the hope that it will be useful,
     152#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     153#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     154#    General Public License for more details.
     155#
     156#  The full text of the GPL is available at:
     157#
     158#                  http://www.gnu.org/licenses/
     159#*****************************************************************************
     160
     161from sage.rings.all import RR, Integer
     162from sage.misc.latex import latex
     163from sage.symbolic.function import BuiltinFunction, is_inexact
     164from sage.symbolic.expression import Expression
     165from sage.structure.coerce import parent
     166import sage.structure.element
     167from sage.calculus.calculus import maxima
     168from sage.libs.mpmath import utils as mpmath_utils
     169from sage.functions.other import sqrt
     170from sage.functions.log import exp
     171from sage.functions.hyperbolic import sinh, cosh
     172from sage.symbolic.constants import pi
     173
     174class Function_Bessel_J(BuiltinFunction):
     175    r"""
     176    The Bessel J Function, denoted by bessel_J(`\alpha`, x) or `J_\alpha(x)`.
     177    As a Taylor series about `x=0` it is equal to:
     178
     179    .. math::
     180
     181        J_\alpha(x) = \sum_{k=0}^\infty \frac{(-1)^k}{k! \Gamma(k+\alpha+1)} (\frac{x}{2})^{2k+\alpha}
     182
     183    For integer orders `\alpha = n` there is an integral representation:
     184
     185    .. math::
     186
     187        J_n(x) = \frac{1}{\pi} \int_0^\pi \cos(n t - x \sin(t)) \; dt
     188
     189    This function also arrises as a special case of the hypergeometric
     190    function `{}_0F_1`:
     191
     192    .. math::
     193
     194        J_\alpha(x) = \frac{x^n}{2^\alpha \Gamma(\alpha + 1)} {}_0F_1(\alpha + 1, -\frac{x^2}{4}).
     195
     196    EXAMPLES::
     197
     198        sage: bessel_J(1, x)
     199        bessel_J(1, x)
     200        sage: bessel_J(1.0, 1.0)
     201        0.440050585744933
     202        sage: n = var('n')
     203        sage: bessel_J(n, x)
     204        bessel_J(n, x)
     205        sage: bessel_J(2, I).n(digits=30)
     206        -0.135747669767038281182852569995
     207
     208    Examples of symbolic manipulation::
     209
     210        sage: a = bessel_J(pi, bessel_J(1, I))
     211        sage: N(a, digits=20)
     212        0.00059023706363796717363 - 0.0026098820470081958110*I
     213
     214        sage: f = bessel_J(2, x)
     215        sage: f.diff(x)
     216        1/2*bessel_J(1, x) - 1/2*bessel_J(3, x)
     217
     218    Currently, integration is not supported (directly) since we cannot yet convert
     219    hypergeometric functions to and from Maxima::
     220
     221        sage: f = bessel_J(2, x)
     222        sage: f.integrate(x)
     223        Traceback (most recent call last):
     224        ...
     225        TypeError: cannot coerce arguments: no canonical coercion from <type 'list'> to Symbolic Ring
     226
     227        sage: m = maxima(bessel_J(2, x))
     228        sage: m.integrate(x)
     229        hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24
     230
     231    Visualization::
     232
     233        sage: plot(bessel_J(1,x), (x,0,5), color='blue')
     234        sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5)) # long time
     235
     236    ALGORITHM:
     237
     238        Numerical evaluation is handled by the mpmath library. Symbolics are
     239        handled by a combination of Maxima and Sage (Ginac/Pynac).
     240    """
     241    def __init__(self):
     242        """
     243        See the docstring for :meth:`Function_Bessel_J`.
     244
     245        EXAMPLES::
     246
     247            sage: sage.functions.bessel.Function_Bessel_J()
     248            bessel_J
     249        """
     250        BuiltinFunction.__init__(self, "bessel_J", nargs=2,
     251                                 conversions=dict(maxima='bessel_j', mathematica='BesselJ'))
     252
     253    def _eval_(self, n, x):
     254        """
     255        EXAMPLES::
     256
     257            sage: a, b = var('a, b')
     258            sage: bessel_J(a, b)
     259            bessel_J(a, b)
     260            sage: bessel_J(1.0, 1.0)
     261            0.440050585744933
     262        """
     263        if not isinstance(n, Expression) and not isinstance(x, Expression) and \
     264               (is_inexact(n) or is_inexact(x)):
     265            coercion_model = sage.structure.element.get_coercion_model()
     266            n, x = coercion_model.canonical_coercion(n, x)
     267            return self._evalf_(n, x, parent(n))
     268
     269        return None # leaves the expression unevaluated
     270
     271    def _evalf_(self, n, x, parent=None):
     272        """
     273        EXAMPLES::
     274
     275            sage: bessel_J(0.0, 1.0)
     276            0.765197686557966
     277            sage: bessel_J(0, 1).n(digits=20)
     278            0.76519768655796655145
     279        """
     280        import mpmath
     281        return mpmath_utils.call(mpmath.besselj, n, x, parent=parent)
     282
     283    def _derivative_(self, n, x, diff_param=None):
     284        """
     285        Return the derivative of the Bessel J function.
     286
     287        EXAMPLES::
     288
     289            sage: f(z) = bessel_J(10, z)
     290            sage: derivative(f, z)
     291            z |--> 1/2*bessel_J(9, z) - 1/2*bessel_J(11, z)
     292        """
     293        return (bessel_J(n-1, x) - bessel_J(n+1, x))/Integer(2)
     294
     295    def _print_latex_(self, n, z):
     296        """
     297        Custom _print_latex_ method.
     298
     299        EXAMPLES::
     300
     301            sage: latex(bessel_J(1, x))
     302            \operatorname{J_{1}}(x)
     303        """
     304        return r"\operatorname{J_{%s}}(%s)" % (latex(n), latex(z))
     305
     306bessel_J = Function_Bessel_J()
     307
     308class Function_Bessel_Y(BuiltinFunction):
     309    r"""
     310    The Bessel Y function.
     311
     312    DEFINITION:
     313
     314    .. math:
     315
     316        bessel\_Y(n, z) = Y_n(z) = ...
     317
     318    The derivative with respect to `z` is:
     319
     320    ..math::
     321
     322        \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\)
     323
     324    EXAMPLES::
     325
     326        sage: bessel_Y(1, x)
     327        bessel_Y(1, x)
     328        sage: bessel_Y(1.0, 1.0)
     329        -0.781212821300289
     330        sage: n = var('n')
     331        sage: bessel_Y(n, x)
     332        bessel_Y(n, x)
     333        sage: bessel_Y(2, I).n()
     334        1.03440456978312 - 0.135747669767038*I
     335
     336    Examples of symbolic manipulation::
     337
     338        sage: a = bessel_Y(pi, bessel_Y(1, I)); a
     339        bessel_Y(pi, bessel_Y(1, I))
     340        sage: N(a, digits=20)
     341        4.2059146571791095708 + 21.307914215321993526*I
     342
     343        sage: f = bessel_Y(2, x)
     344        sage: f.diff(x)
     345        1/2*bessel_Y(1, x) - 1/2*bessel_Y(3, x)
     346
     347    High preceision and complex valued inputs (fixes Trac issue #4230)::
     348
     349        sage: bessel_Y(0, 1).n(128)
     350        0.088256964215676957982926766023515162828
     351        sage: bessel_Y(0, RealField(200)(1))
     352        0.088256964215676957982926766023515162827817523090675546711044
     353        sage: bessel_Y(0, ComplexField(200)(0.5+I))
     354        0.077763160184438051408593468823822434235010300228009867784073 + 1.0142336049916069152644677682828326441579314239591288411739*I
     355
     356    Visualization::
     357
     358        sage: plot(bessel_Y(1,x), (x,0,5), color='blue')
     359        sage: complex_plot(bessel_Y(1, x), (-5, 5), (-5, 5)) # long time
     360
     361    ALGORITHM:
     362
     363        Numerical evaluation is handled by the mpmath library. Symbolics are
     364        handled by a combination of Maxima and Sage (Ginac/Pynac).
     365    """
     366    def __init__(self):
     367        """
     368        See the docstring for :meth:`Function_Bessel_Y`.
     369
     370        EXAMPLES::
     371
     372            sage: sage.functions.bessel.Function_Bessel_Y()(0, x)
     373            bessel_Y(0, x)
     374        """
     375        BuiltinFunction.__init__(self, "bessel_Y", nargs=2,
     376                                 conversions=dict(maxima='bessel_y', mathematica='BesselY'))
     377
     378    def _eval_(self, n, x):
     379        """
     380        EXAMPLES::
     381
     382            sage: a,b = var('a, b')
     383            sage: bessel_Y(a, b)
     384            bessel_Y(a, b)
     385            sage: bessel_Y(0, 1).n(128)
     386            0.088256964215676957982926766023515162828
     387        """
     388        if not isinstance(n, Expression) and not isinstance(x, Expression) and \
     389               (is_inexact(n) or is_inexact(x)):
     390            coercion_model = sage.structure.element.get_coercion_model()
     391            n, x = coercion_model.canonical_coercion(n, x)
     392            return self._evalf_(n, x, parent(n))
     393
     394        return None # leaves the expression unevaluated
     395
     396    def _evalf_(self, n, x, parent=None):
     397        """
     398        EXAMPLES::
     399
     400            sage: bessel_Y(1.0+2*I, 3.0+4*I)
     401            0.699410324467538 + 0.228917940896421*I
     402            sage: bessel_Y(0, 1).n(256)
     403            0.08825696421567695798292676602351516282781752309067554671104384761199978932351
     404        """
     405        import mpmath
     406        return mpmath_utils.call(mpmath.bessely, n, x, parent=parent)
     407
     408    def _derivative_(self, n, x, diff_param=None):
     409        """
     410        Return the derivative of the Bessel Y function.
     411
     412        EXAMPLES::
     413
     414            sage: f(x) = bessel_Y(10, x)
     415            sage: derivative(f, x)
     416            x |--> 1/2*bessel_Y(9, x) - 1/2*bessel_Y(11, x)
     417        """
     418        return (bessel_Y(n-1, x) - bessel_Y(n+1, x))/Integer(2)
     419
     420    def _print_latex_(self, n, z):
     421        """
     422        Custom _print_latex_ method.
     423
     424        EXAMPLES::
     425
     426            sage: latex(bessel_Y(1, x))
     427            \operatorname{Y_{1}}(x)
     428        """
     429        return r"\operatorname{Y_{%s}}(%s)" % (latex(n), latex(z))
     430
     431bessel_Y = Function_Bessel_Y()
     432
     433class Function_Bessel_I(BuiltinFunction):
     434    r"""
     435    The Bessel I function, or the Modified Bessel Function of the First Kind.
     436
     437    DEFINITION:
     438
     439    .. math::
     440
     441        \operatorname{bessel\_I}(\alpha, x) = I_\alpha(x) = i^{-\alpha} J_\alpha(ix)
     442
     443    EXAMPLES::
     444
     445        sage: bessel_I(1, x)
     446        bessel_I(1, x)
     447        sage: bessel_I(1.0, 1.0)
     448        0.565159103992485
     449        sage: n = var('n')
     450        sage: bessel_I(n, x)
     451        bessel_I(n, x)
     452        sage: bessel_I(2, I).n()
     453        -0.114903484931900
     454
     455    Examples of symbolic manipulation::
     456
     457        sage: a = bessel_I(pi, bessel_I(1, I))
     458        sage: N(a, digits=20)
     459        0.00026073272117205890528 - 0.0011528954889080572266*I
     460
     461        sage: f = bessel_I(2, x)
     462        sage: f.diff(x)
     463        1/2*bessel_I(1, x) + 1/2*bessel_I(3, x)
     464
     465    Special identities that bessel_I satisfies::
     466
     467        sage: bessel_I(1/2, x)
     468        sqrt(1/(pi*x))*sqrt(2)*sinh(x)
     469        sage: eq = bessel_I(1/2, x) == bessel_I(0.5, x)
     470        sage: eq.test_relation()
     471        True
     472        sage: bessel_I(-1/2, x)
     473        sqrt(1/(pi*x))*sqrt(2)*cosh(x)
     474        sage: eq = bessel_I(-1/2, x) == bessel_I(-0.5, x)
     475        sage: eq.test_relation()
     476        True
     477
     478    High preceision and complex valued inputs::
     479
     480        sage: bessel_I(0, 1).n(128)
     481        1.2660658777520083355982446252147175376
     482        sage: bessel_I(0, RealField(200)(1))
     483        1.2660658777520083355982446252147175376076703113549622068081
     484        sage: bessel_I(0, ComplexField(200)(0.5+I))
     485        0.80644357583493619472428518415019222845373366024179916785502 + 0.22686958987911161141397453401487525043310874687430711021434*I
     486
     487    Visualization::
     488
     489        sage: plot(bessel_I(1,x), (x,0,5), color='blue')
     490        sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5)) # long time
     491
     492    ALGORITHM:
     493
     494        Numerical evaluation is handled by the mpmath library. Symbolics are
     495        handled by a combination of Maxima and Sage (Ginac/Pynac).
     496    """
     497    def __init__(self):
     498        """
     499        See the docstring for :meth:`Function_Bessel_I`.
     500
     501        EXAMPLES::
     502
     503            sage: bessel_I(1,x)
     504            bessel_I(1, x)
     505        """
     506        BuiltinFunction.__init__(self, "bessel_I", nargs=2,
     507                                 conversions=dict(maxima='bessel_i', mathematica='BesselI'))
     508
     509    def _eval_(self, n, x):
     510        """
     511        EXAMPLES::
     512
     513            sage: y=var('y')
     514            sage: bessel_I(y,x)
     515            bessel_I(y, x)
     516            sage: bessel_I(0.0, 1.0)
     517            1.26606587775201
     518            sage: bessel_I(1/2, 1)
     519            sqrt(2)*sinh(1)/sqrt(pi)
     520            sage: bessel_I(-1/2, pi)
     521            sqrt(2)*cosh(pi)/pi
     522        """
     523        if not isinstance(n, Expression) and not isinstance(x, Expression) and \
     524               (is_inexact(n) or is_inexact(x)):
     525            coercion_model = sage.structure.element.get_coercion_model()
     526            n, x = coercion_model.canonical_coercion(n, x)
     527            return self._evalf_(n, x, parent(n))
     528
     529        # special identities
     530        if n == Integer(1)/Integer(2):
     531            return sqrt(2/(pi*x))*sinh(x)
     532        elif n == -Integer(1)/Integer(2):
     533            return sqrt(2/(pi*x))*cosh(x)
     534
     535        return None # leaves the expression unevaluated
     536
     537    def _evalf_(self, n, x, parent=None):
     538        """
     539        EXAMPLES::
     540
     541            sage: bessel_I(1,3).n(digits=20)
     542            3.9533702174026093965
     543        """
     544        import mpmath
     545        return mpmath_utils.call(mpmath.besseli, n, x, parent=parent)
     546
     547    def _derivative_(self, n, x, diff_param=None):
     548        """
     549        Return the derivative of the Bessel I function `I_n(x)` with repect
     550        to `x`.
     551
     552        EXAMPLES::
     553
     554            sage: f(z) = bessel_I(10, x)
     555            sage: derivative(f, x)
     556            z |--> 1/2*bessel_I(9, x) + 1/2*bessel_I(11, x)
     557        """
     558        return (bessel_I(n-1, x) + bessel_I(n+1, x))/Integer(2)
     559
     560    def _print_latex_(self, n, z):
     561        """
     562        Custom _print_latex_ method.
     563
     564        EXAMPLES::
     565
     566            sage: latex(bessel_I(1, x))
     567            \operatorname{I_{1}}(x)
     568        """
     569        return r"\operatorname{I_{%s}}(%s)" % (latex(n), latex(z))
     570
     571bessel_I = Function_Bessel_I()
     572
     573class Function_Bessel_K(BuiltinFunction):
     574    r"""
     575    The Bessel K function.
     576
     577    DEFINITION:
     578
     579    .. math:
     580
     581        \operatorname{bessel\_K}(\alpha, x) = K_\alpha(x) = ...
     582
     583    EXAMPLES::
     584
     585        sage: bessel_K(1, x)
     586        bessel_K(1, x)
     587        sage: bessel_K(1.0, 1.0)
     588        0.601907230197235
     589        sage: n = var('n')
     590        sage: bessel_K(n, x)
     591        bessel_K(n, x)
     592        sage: bessel_K(2, I).n()
     593        -2.59288617549120 + 0.180489972066962*I
     594
     595    Examples of symbolic manipulation::
     596
     597        sage: a = bessel_K(pi, bessel_K(1, I)); a
     598        bessel_K(pi, bessel_K(1, I))
     599        sage: N(a, digits=20)
     600        3.8507583115005220157 + 0.068528298579883425792*I
     601
     602        sage: f = bessel_K(2, x)
     603        sage: f.diff(x)
     604        1/2*bessel_K(1, x) + 1/2*bessel_K(3, x)
     605
     606        sage: bessel_K(1/2, x)
     607        bessel_K(1/2, x)
     608        sage: bessel_K(1/2, -1)
     609        bessel_K(1/2, -1)
     610        sage: bessel_K(1/2, 1)
     611        sqrt(pi)*sqrt(1/2)*e^(-1)
     612
     613    High preceision and complex valued inputs::
     614
     615        sage: bessel_K(0, 1).n(128)
     616        0.42102443824070833333562737921260903614
     617        sage: bessel_K(0, RealField(200)(1))
     618        0.42102443824070833333562737921260903613621974822666047229897
     619        sage: bessel_K(0, ComplexField(200)(0.5+I))
     620        0.058365979093103864080375311643360048144715516692187818271179 - 0.67645499731334483535184142196073004335768129348518210260256*I
     621
     622    Visualization::
     623
     624        sage: plot(bessel_K(1,x), (x,0,5), color='blue')
     625        sage: complex_plot(bessel_K(1, x), (-5, 5), (-5, 5)) # long time
     626
     627    ALGORITHM:
     628
     629        Numerical evaluation is handled by the mpmath library. Symbolics are
     630        handled by a combination of Maxima and Sage (Ginac/Pynac).
     631
     632    TESTS:
     633
     634    Verfify that Trac #3426 is fixed:
     635
     636    The Bessel K function can be evaluated numerically at complex orders::
     637
     638        sage: bessel_K(10 * I, 10).n()
     639        9.82415743819925e-8
     640
     641    For a fixed imaginary order and increasin, real, second component the
     642    value of Bessel K is exponentially decaying::
     643
     644        sage: for x in [10, 20, 50, 100, 200]: print bessel_K(5*I, x).n()
     645        5.27812176514912e-6
     646        3.11005908421801e-10
     647        2.66182488515423e-23 - 8.59622057747552e-58*I
     648        4.11189776828337e-45 - 1.01494840019482e-80*I
     649        1.15159692553603e-88 - 6.75787862113718e-125*I
     650    """
     651    def __init__(self):
     652        """
     653        See the docstring for :meth:`Function_Bessel_K`.
     654
     655        EXAMPLES::
     656
     657            sage: sage.functions.bessel.Function_Bessel_K()
     658            bessel_K
     659        """
     660        BuiltinFunction.__init__(self, "bessel_K", nargs=2,
     661                                 conversions=dict(maxima='bessel_k', mathematica='BesselK'))
     662
     663    def _eval_(self, n, x):
     664        """
     665        EXAMPLES::
     666
     667            sage: bessel_K(1,0)
     668            bessel_K(1, 0)
     669            sage: bessel_K(1.0, 0.0)
     670            +infinity
     671            sage: bessel_K(-1, 1).n(128)
     672            0.60190723019723457473754000153561733926
     673        """
     674        if not isinstance(n, Expression) and not isinstance(x, Expression) and \
     675               (is_inexact(n) or is_inexact(x)):
     676            coercion_model = sage.structure.element.get_coercion_model()
     677            n, x = coercion_model.canonical_coercion(n, x)
     678            return self._evalf_(n, x, parent(n))
     679
     680        # special identity
     681        if n == Integer(1)/Integer(2) and x > 0:
     682            return sqrt(pi/2)*exp(-x)*x**(Integer(1)/Integer(2))
     683
     684        return None # leaves the expression unevaluated
     685
     686    def _evalf_(self, n, x, parent=None):
     687        """
     688        EXAMPLES::
     689
     690            sage: bessel_K(0.0, 1.0)
     691            0.421024438240708
     692            sage: bessel_K(0, RealField(128)(1))
     693            0.42102443824070833333562737921260903614
     694        """
     695        import mpmath
     696        return mpmath_utils.call(mpmath.besselk, n, x, parent=parent)
     697
     698    def _derivative_(self, n, x, diff_param=None):
     699        """
     700        Return the derivative of the Bessel K function.
     701
     702        EXAMPLES::
     703
     704            sage: f(x) = bessel_K(10, x)
     705            sage: derivative(f, x)
     706            x |--> 1/2*bessel_K(9, x) + 1/2*bessel_K(11, x)
     707        """
     708        return (bessel_K(n-1, x) + bessel_K(n+1, x))/Integer(2)
     709
     710    def _print_latex_(self, n, z):
     711        """
     712        Custom _print_latex_ method.
     713
     714        EXAMPLES::
     715
     716            sage: latex(bessel_K(1, x))
     717            \operatorname{K_{1}}(x)
     718        """
     719        return r"\operatorname{K_{%s}}(%s)" % (latex(n), latex(z))
     720
     721bessel_K = Function_Bessel_K()
     722
     723
     724# dictionary used in Bessel
     725bessel_type_dict = {'I': bessel_I, 'J': bessel_J, 'K': bessel_K, 'Y': bessel_Y}
     726
     727def Bessel(*args, **kwds):
     728    """
     729    A factory function that produces symbolic I, J, K, and Y Bessel functions. There
     730    are several ways to call this function:
     731
     732        - ``Bessel(order, type)``
     733        - ``Bessel(order)`` -- type defaults to 'J'
     734        - ``Bessel(order, typ=T)``
     735        - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter function
     736        - ``Bessel()`` -- order is unspecified, type is 'J'
     737
     738    where `order` can be any integer and T must be one of the strings 'I', 'J',
     739    'K', or 'Y'.
     740
     741    See the EXAMPLES below.
     742
     743    EXAMPLES:
     744
     745    Construction of Bessel functions with various orders and types::
     746
     747        sage: Bessel()
     748        bessel_J
     749        sage: Bessel(1)(x)
     750        bessel_J(1, x)
     751        sage: Bessel(1, 'Y')(x)
     752        bessel_Y(1, x)
     753        sage: Bessel(-2, 'Y')(x)
     754        bessel_Y(-2, x)
     755        sage: Bessel(typ='K')
     756        bessel_K
     757        sage: Bessel(0, typ='I')(x)
     758        bessel_I(0, x)
     759
     760    Evaluation::
     761
     762        sage: f = Bessel(1)
     763        sage: f(3.0)
     764        0.339058958525936
     765        sage: f(3)
     766        bessel_J(1, 3)
     767        sage: f(3).n(digits=50)
     768        0.33905895852593645892551459720647889697308041819801
     769
     770        sage: g = Bessel(typ='J')
     771        sage: g(1,3)
     772        bessel_J(1, 3)
     773        sage: g(2, 3+I).n()
     774        0.634160370148554 + 0.0253384000032695*I
     775        sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15
     776        True
     777
     778    Symbolic calculus::
     779
     780        sage: f(x) = Bessel(0, 'J')(x)
     781        sage: derivative(f, x)
     782        x |--> 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x)
     783        sage: derivative(f, x, x)
     784        x |--> 1/4*bessel_J(-2, x) - 1/2*bessel_J(0, x) + 1/4*bessel_J(2, x)
     785
     786    Verify that `J_0` satisfies Bessel's differential equation numerically
     787    using the `test_relation()` method::
     788
     789        sage: y = bessel_J(0, x)
     790        sage: diffeq = x^2*derivative(y,x,x) + x*derivative(y,x) + x^2*y == 0
     791        sage: diffeq.test_relation(proof=False)
     792        True
     793
     794    Conversion to other systems::
     795
     796        sage: x,y = var('x,y')
     797        sage: f = maxima(Bessel(typ='K')(x,y))
     798        sage: f.derivative('x')
     799        %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)
     800        sage: f.derivative('y')
     801        -(bessel_k(x+1,y)+bessel_k(x-1,y))/2
     802
     803    Compute the particular solution to Bessel's Differential Equation that
     804    satisfies `y(1) = 1` and `y'(1) = 1`, then verify the initial conditions
     805    and plot it::
     806
     807        sage: y = function('y', x)
     808        sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0
     809        sage: f = desolve(diffeq, y, [1, 1, 1]); f
     810        (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) - (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1))
     811
     812        sage: f.subs(x=1).n() # numerical verification
     813        1.00000000000000
     814        sage: fp = f.diff(x)
     815        sage: fp.subs(x=1).n()
     816        1.00000000000000
     817
     818        sage: f.subs(x=1).simplify_full() # symbolic verification
     819        1
     820        sage: fp = f.diff(x)
     821        sage: fp.subs(x=1).simplify_full()
     822        1
     823
     824        sage: plot(f, (x,0,5))
     825
     826    Plotting::
     827
     828        sage: f(x) = Bessel(0)(x); f
     829        x |--> bessel_J(0, x)
     830        sage: plot(f, (x, 1, 10))
     831
     832        sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10)
     833
     834        sage: G = Graphics()
     835        sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ])
     836        sage: show(G)
     837
     838    A recreation of Abramowitz and Stegun Figure 9.1::
     839
     840        sage: G  = plot(Bessel(0, 'J'), 0, 15, color='black')
     841        sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black')
     842        sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted')
     843        sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted')
     844        sage: show(G, ymin=-1, ymax=1)
     845
     846    """
     847    # Determine the order and type of function from the arguments and keywords.
     848    # These are recored in local variables: _type, _order, _system, _nargs.
     849    _type = None
     850    if len(args) == 0: # no order specified
     851        _order = None
     852        _nargs = 2
     853    elif len(args) == 1: # order is specified
     854        _order = args[0]
     855        _nargs = 1
     856    elif len(args) == 2: # both order and type are positional arguments
     857        _order = args[0]
     858        _type = args[1]
     859        _nargs = 1
     860    else:
     861        raise TypeError("at most two positional arguments may be specified, "
     862                       +"see the docstring for Bessel")
     863    # check for type inconsistency
     864    if _type is not None and 'typ' in kwds and _type != kwds['typ']:
     865        raise ValueError("inconsistent types given")
     866    # record the function type
     867    if _type is None:
     868        if 'typ' in kwds:
     869            _type = kwds['typ']
     870        else:
     871            _type = 'J'
     872    if not (_type in ['I', 'J', 'K', 'Y']):
     873        raise ValueError("type must be one of I, J, K, Y")
     874    # record the numerical evaluation system
     875    if 'algorithm' in kwds:
     876        _system = kwds['algorithm']
     877    else:
     878        _system = 'mpmath'
     879
     880    # return the function
     881    # TODO remove assertions
     882    assert _type in ['I', 'J', 'K', 'Y']
     883    assert _order is None or _order in RR
     884    assert _nargs == 1 or _nargs == 2
     885    assert _system == 'mpmath'
     886
     887    # TODO what to do with _system?
     888    _f = bessel_type_dict[_type]
     889    if _nargs == 1:
     890        return lambda x: _f(_order, x)
     891    else:
     892        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.