Ticket #12455: trac_12455-newstyle-airy-rebase.patch

File trac_12455-newstyle-airy-rebase.patch, 31.7 KB (added by kcrisman, 6 years ago)
  • new file sage/functions/airy.py

    # HG changeset patch
    # User D. S. McNeil <dsm054@gmail.com>
    # Date 1338267871 25200
    # Node ID c062930f0159340f9d200948bfbacc0aa4def630
    # Parent  bd52cf96d40a630c7c313990861141658882ee96
    Trac 12455 - Convert Airy functions to symbolic
    
    diff --git a/sage/functions/airy.py b/sage/functions/airy.py
    new file mode 100644
    - +  
     1#*****************************************************************************
     2#       Copyright (C) 2010 Oscar Gerardo Lazo Arjona algebraicamente@gmail.com
     3#
     4#  Distributed under the terms of the GNU General Public License (GPL)
     5#
     6#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     7#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     8#    General Public License for more details.
     9#
     10#  The full text of the GPL is available at:
     11#
     12#                  http://www.gnu.org/licenses/
     13#*****************************************************************************
     14
     15from sage.symbolic.function import BuiltinFunction
     16from sage.symbolic.expression import Expression
     17from sage.symbolic.function import is_inexact
     18from sage.structure.coerce import parent as sage_structure_coerce_parent
     19
     20from sage.functions.other import gamma
     21from sage.rings.integer_ring import ZZ
     22from sage.rings.real_double import RDF
     23from sage.functions.special import meval
     24from sage.calculus.var import var
     25from sage.calculus.functional import derivative
     26
     27# from sage.rings.real_mpfr import RR
     28from sage.rings.rational import Rational as R
     29
     30class FunctionAiryAiGeneral(BuiltinFunction):
     31    def __init__(self):
     32        r"""
     33        The generalized derivative of the Airy Ai function
     34
     35        INPUT:
     36       
     37        - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`.
     38          For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Ai}^{(n)}(z)`,
     39          and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
     40         
     41        .. math ::
     42         
     43            f_0(z) = \operatorname{Ai}(z)
     44         
     45            f_n(z) = \int_0^z f_{n-1}(t) dt.
     46
     47        - ``x``.- The argument of the function.
     48
     49        - ``hold``.- Whether or not to stop from returning higher derivatives in terms of
     50          `\operatorname{Ai}(x)` and `\operatorname{Ai}'(x)`.
     51       
     52        EXAMPLES::
     53
     54            sage: x,n=var('x n')
     55            sage: airy_ai(-2,x)
     56            airy_ai(-2, x)
     57            sage: derivative(airy_ai(-2,x),x)
     58            airy_ai(-1, x)
     59            sage: airy_ai(n,x)
     60            airy_ai(n, x)
     61            sage: derivative(airy_ai(n,x),x)
     62            airy_ai(n + 1, x)
     63            sage: airy_ai(2,x,True)
     64            airy_ai(2, x)
     65            sage: derivative(airy_ai(2,x,hold_derivative=True),x)
     66            airy_ai(3, x)
     67        """
     68
     69        BuiltinFunction.__init__(self, "airy_ai", nargs=2,
     70                latex_name=r"\operatorname{Ai}")
     71
     72    def _derivative_(self, alpha, *args, **kwds):
     73        """
     74        EXAMPLES::
     75
     76            sage: x,n=var('x n')
     77            sage: derivative(airy_ai(n,x),x) # indirect doctest
     78            airy_ai(n + 1, x)
     79        """
     80
     81        x=args[0]
     82        return airy_ai_general(alpha+1,x)
     83
     84    def _eval_(self, alpha, *args):
     85        """
     86        EXAMPLES::
     87
     88            sage: x,n=var('x n')
     89            sage: airy_ai(-2,1.0) # indirect doctest
     90            0.136645379421096
     91            sage: airy_ai(n,1.0) # indirect doctest
     92            airy_ai(n, 1.00000000000000)
     93        """
     94
     95        x=args[0]
     96        if not isinstance(x,Expression) and not isinstance(alpha,Expression):
     97            if is_inexact(x):
     98                return self._evalf_(alpha,x)
     99        else:
     100            return None
     101
     102    def _evalf_(self, alpha, x, **kwargs):
     103        """
     104        EXAMPLES::
     105
     106            sage: airy_ai(-2,1.0) # indirect doctest
     107            0.136645379421096           
     108        """
     109        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     110        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     111        prec = parent.prec() if hasattr(parent, 'prec') else 53
     112        if algorithm == 'mpmath':
     113            import mpmath
     114            from sage.libs.mpmath import utils as mpmath_utils
     115            return mpmath_utils.call(mpmath.airyai, x, derivative=alpha, parent=parent)
     116        elif algorithm == 'maxima':
     117            raise NotImplementedError("general case not available in maxima")
     118        else:
     119            raise ValueError("unknown algorithm")
     120
     121class FunctionAiryAiSimple(BuiltinFunction):
     122    def __init__(self):
     123        """
     124        The class for the Airy Ai function.
     125
     126        EXAMPLES::
     127
     128            sage: from sage.functions.airy import FunctionAiryAiSimple
     129            sage: airy_ai_simple= FunctionAiryAiSimple()
     130            sage: f=airy_ai_simple(x); f
     131            airy_ai(x)
     132        """
     133
     134        BuiltinFunction.__init__(self, "airy_ai", latex_name=r'\operatorname{Ai}',
     135                conversions=dict(mathematica='AiryAi', maxima='airy_ai'))
     136
     137    def _derivative_(self, x, diff_param=None):
     138        """
     139        EXAMPLES::
     140       
     141            sage: derivative(airy_ai(x),x) # indirect doctest
     142            airy_ai_prime(x)
     143        """
     144        return airy_ai_prime(x)
     145
     146    def _eval_(self, x):
     147        """
     148        EXAMPLES::
     149       
     150            sage: airy_ai(0) # indirect doctest
     151            1/3*3^(1/3)/gamma(2/3)
     152            sage: airy_ai(0.0) # indirect doctest
     153            0.355028053887817
     154            sage: airy_ai(I) # indirect doctest
     155            airy_ai(I)
     156            sage: airy_ai(1.0*I) # indirect doctest
     157            0.331493305432141 - 0.317449858968444*I
     158        """
     159
     160        if not isinstance(x,Expression):
     161            if is_inexact(x):
     162                return self._evalf_(x)
     163            elif x==0:
     164                r=ZZ(2)/3
     165                return 1/(3**(r)*gamma(r))
     166        else:
     167            return None
     168
     169    def _evalf_(self, x, **kwargs):
     170        """
     171        EXAMPLES::
     172       
     173            sage: airy_ai(0.0) # indirect doctest
     174            0.355028053887817
     175            sage: airy_ai(1.0*I) # indirect doctest
     176            0.331493305432141 - 0.317449858968444*I
     177
     178        We can use several methods for numerical evaluation::
     179
     180            sage: airy_ai(3).n(algorithm='maxima')
     181            0.00659113935746
     182            sage: airy_ai(3).n(algorithm='mpmath')
     183            0.00659113935746072
     184            sage: airy_ai(3).n(algorithm='mpmath',prec=100)
     185            0.0065911393574607191442574484080
     186           
     187        TESTS::
     188
     189            sage: airy_ai(3).n(algorithm='maxima', prec=70)
     190            Traceback (most recent call last):
     191            ...
     192            ValueError: for the maxima algorithm the precision must be 53
     193       
     194        """
     195        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     196        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     197        prec = parent.prec() if hasattr(parent, 'prec') else 53
     198        if algorithm == 'mpmath':
     199            import mpmath
     200            from sage.libs.mpmath import utils as mpmath_utils
     201            return mpmath_utils.call(mpmath.airyai, x, parent=parent)
     202        elif algorithm == 'maxima':
     203            if prec != 53:
     204                raise ValueError, "for the maxima algorithm the precision must be 53"
     205            return RDF(meval("airy_ai(%s)" % RDF(x)))
     206        else:
     207            raise ValueError("unknown algorithm")
     208
     209class FunctionAiryAiPrime(BuiltinFunction):
     210    def __init__(self):
     211        """
     212        The derivative of the Airy Ai function; see :func:`airy_ai`
     213        for the full documentation.
     214
     215        EXAMPLES::
     216
     217            sage: x,n=var('x n')
     218            sage: airy_ai_prime(x)
     219            airy_ai_prime(x)
     220            sage: airy_ai_prime(0)
     221            -1/3*3^(2/3)/gamma(1/3)
     222        """
     223
     224        BuiltinFunction.__init__(self, "airy_ai_prime",
     225                latex_name=r"\operatorname{Ai}'",
     226                conversions=dict(mathematica='AiryAiPrime', maxima='airy_dai'))
     227
     228    def _derivative_(self, x, diff_param=None):
     229        """
     230        EXAMPLES::
     231
     232            sage: derivative(airy_ai_prime(x),x) # indirect doctest
     233            x*airy_ai(x)
     234        """
     235        return x*airy_ai_simple(x)
     236
     237    def _eval_(self, x):
     238        """
     239        EXAMPLES::
     240
     241            sage: airy_ai(1,0) # indirect doctest
     242            -1/3*3^(2/3)/gamma(1/3)
     243            sage: airy_ai(1,0.0) # indirect doctest
     244            -0.258819403792807
     245        """
     246        if not isinstance(x,Expression):
     247            if is_inexact(x):
     248                return self._evalf_(x)
     249            elif x==0:
     250                r=ZZ(1)/3
     251                return -1/(3**(r)*gamma(r))
     252        else:
     253            return None
     254
     255    def _evalf_(self, x, **kwargs):
     256        """
     257        EXAMPLES::
     258       
     259            sage: airy_ai(1,0.0) # indirect doctest
     260            -0.258819403792807
     261
     262        We can use several methods for numerical evaluation::
     263
     264            sage: airy_ai(1,4).n(algorithm='maxima')
     265            -0.0019586409502
     266            sage: airy_ai(1,4).n(algorithm='mpmath')
     267            -0.00195864095020418
     268            sage: airy_ai(1,4).n(algorithm='mpmath', prec=100)
     269            -0.0019586409502041789001381409184
     270
     271        TESTS::
     272
     273            sage: airy_ai(1, 4).n(algorithm='maxima', prec=70)
     274            Traceback (most recent call last):
     275            ...
     276            ValueError: for the maxima algorithm the precision must be 53
     277
     278
     279        """
     280        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     281        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     282        prec = parent.prec() if hasattr(parent, 'prec') else 53
     283        if algorithm == 'mpmath':
     284            import mpmath
     285            from sage.libs.mpmath import utils as mpmath_utils
     286            return mpmath_utils.call(mpmath.airyai, x, derivative=1, parent=parent)
     287        elif algorithm == 'maxima':
     288            if prec != 53:
     289                raise ValueError, "for the maxima algorithm the precision must be 53"
     290            return RDF(meval("airy_dai(%s)" % RDF(x)))
     291        else:
     292            raise ValueError("unknown algorithm")
     293
     294
     295airy_ai_general=FunctionAiryAiGeneral()
     296airy_ai_simple= FunctionAiryAiSimple()
     297airy_ai_prime=  FunctionAiryAiPrime()
     298
     299
     300def airy_ai(alpha,x=None, hold_derivative=False, *args, **kwds):
     301    r"""
     302    The Airy Ai function `\operatorname{Ai}(x)` is one of the two
     303    linearly independent solutions to the Airy differental equation
     304    `f''(z) +f(z)x=0`, defined by the initial conditions:
     305
     306    .. math ::
     307        \operatorname{Ai}(0)=\frac{1}{2^{2/3} \Gamma(\frac{2}{3})},
     308
     309        \operatorname{Ai}'(0)=-\frac{1}{2^{1/3} \Gamma(\frac{1}{3})}.
     310
     311    Another way to define the Airy Ai function is:
     312
     313    .. math::
     314        \operatorname{Ai}(x)=\frac{1}{\pi}\int_0^\infty
     315        \cos\left(\frac{1}{3}t^3+xt\right) dt.
     316
     317    INPUT:
     318   
     319    - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`.
     320      For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Ai}^{(n)}(z)`,
     321      and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
     322     
     323    .. math ::
     324     
     325        f_0(z) = \operatorname{Ai}(z)
     326     
     327        f_n(z) = \int_0^z f_{n-1}(t) dt.
     328
     329    - ``x``.- The argument of the function.
     330
     331    - ``hold_derivative``.- Whether or not to stop from returning higher derivatives
     332       in terms of `\operatorname{Ai}(x)` and `\operatorname{Ai}'(x)`.
     333   
     334
     335    EXAMPLES::
     336
     337        sage: n,x=var('n x')
     338        sage: airy_ai(x)
     339        airy_ai(x)
     340
     341    It can return derivatives or integrals::
     342   
     343        sage: airy_ai(1,x)
     344        airy_ai_prime(x)
     345        sage: airy_ai(2,x)
     346        x*airy_ai(x)
     347        sage: airy_ai(2,x,True)
     348        airy_ai(2, x)
     349        sage: airy_ai(-2,x)
     350        airy_ai(-2, x)
     351        sage: airy_ai(n, x)
     352        airy_ai(n, x)
     353
     354    It can be evaluated symbolically or numerically for real or complex values::
     355
     356        sage: airy_ai(0)
     357        1/3*3^(1/3)/gamma(2/3)
     358        sage: airy_ai(0.0)
     359        0.355028053887817
     360        sage: airy_ai(I)
     361        airy_ai(I)
     362        sage: airy_ai(1.0*I)
     363        0.331493305432141 - 0.317449858968444*I
     364
     365    The functions can be evaluated numerically using either mpmath (the default)
     366    or maxima, where mpmath can compute the values to arbitrary precision::
     367
     368        sage: airy_ai(2).n(prec=100)
     369        0.034924130423274379135322080792
     370        sage: airy_ai(2).n(algorithm='mpmath',prec=100)
     371        0.034924130423274379135322080792
     372        sage: airy_ai(2).n(algorithm='maxima')
     373        0.0349241304233
     374
     375    And the derivatives can be evaluated::
     376
     377        sage: airy_ai(1,0)
     378        -1/3*3^(2/3)/gamma(1/3)
     379        sage: airy_ai(1,0.0)
     380        -0.258819403792807
     381
     382    Plots::
     383
     384        sage: plot(airy_ai(x),(x,-10,5))+plot(airy_ai_prime(x),(x,-10,5),color='red')
     385       
     386    **References**
     387     
     388    - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10"
     389
     390    - http://en.wikipedia.org/wiki/Airy_function
     391    """
     392
     393    #We catch the case with no alpha
     394    if x==None:
     395        x=alpha
     396        return airy_ai_simple(x, **kwds)
     397    #We raise an error if there are too many arguments
     398    if len(args) > 0:
     399        raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" %
     400                        (len(args)+3))
     401
     402    #We take care of all other cases.
     403    if hold_derivative:
     404        return airy_ai_general(alpha,x,**kwds)
     405    elif alpha==0:
     406        return airy_ai_simple(x, **kwds)
     407    elif alpha==1:
     408        return airy_ai_prime(x, **kwds)
     409    elif alpha>1:
     410        #we use a different variable here because if x is a
     411        #particular value, we would be differentiating a constant
     412        #which would return 0. What we want is the value of
     413        #the derivative at the value and not the derivative of
     414        #a particular value of the function.
     415        v=var('v')
     416        return derivative(airy_ai_simple(v,**kwds),v,alpha).subs(v=x)
     417    else:
     418        return airy_ai_general(alpha,x,**kwds)
     419
     420########################################################################
     421########################################################################
     422
     423class FunctionAiryBiGeneral(BuiltinFunction):
     424    def __init__(self):
     425        r"""
     426        The generalized derivative of the Airy Bi function
     427
     428        INPUT:
     429       
     430        - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to
     431          `z`. For `\alpha = n = 1,2,3,\ldots` this gives the derivative
     432          `\operatorname{Bi}^{(n)}(z)`, and for `\alpha = -n = -1,-2,-3,\ldots` this gives
     433          the `n`-fold iterated integral
     434         
     435        .. math ::
     436         
     437            f_0(z) = \operatorname{Bi}(z)
     438         
     439            f_n(z) = \int_0^z f_{n-1}(t) dt.
     440
     441        - ``x``.- The argument of the function.
     442
     443        - ``hold``.- Whether or not to stop from returning higher derivatives in terms of
     444          `\operatorname{Bi}(x)` and `\operatorname{Bi}'(x)`.
     445       
     446        EXAMPLES::
     447
     448            sage: x,n=var('x n')
     449            sage: airy_bi(-2,x)
     450            airy_bi(-2, x)
     451            sage: derivative(airy_bi(-2,x),x)
     452            airy_bi(-1, x)
     453            sage: airy_bi(n,x)
     454            airy_bi(n, x)
     455            sage: derivative(airy_bi(n,x),x)
     456            airy_bi(n + 1, x)
     457            sage: airy_bi(2,x,True)
     458            airy_bi(2, x)
     459            sage: derivative(airy_bi(2,x,hold_derivative=True),x)
     460            airy_bi(3, x)
     461        """
     462
     463        BuiltinFunction.__init__(self, "airy_bi", nargs=2,
     464            latex_name=r"\operatorname{Bi}")
     465
     466    def _derivative_(self, alpha, *args, **kwds):
     467        """
     468        EXAMPLES::
     469
     470            sage: x,n=var('x n')
     471            sage: derivative(airy_bi(n,x),x) # indirect doctest
     472            airy_bi(n + 1, x)
     473        """
     474
     475        x=args[0]
     476        return airy_bi_general(alpha+1,x)
     477
     478    def _eval_(self, alpha, *args):
     479        """
     480        EXAMPLES::
     481
     482            sage: x,n=var('x n')
     483            sage: airy_bi(-2,1.0) # indirect doctest
     484            0.388621540699059
     485            sage: airy_bi(n,1.0) # indirect doctest
     486            airy_bi(n, 1.00000000000000)
     487        """
     488
     489       
     490        x=args[0]
     491        if not isinstance(x,Expression) and not isinstance(alpha,Expression):
     492            if is_inexact(x):
     493                return self._evalf_(alpha,x)
     494        else:
     495            return None
     496
     497    def _evalf_(self, alpha, x, **kwargs):
     498        """
     499        EXAMPLES::
     500
     501            sage: airy_bi(-2,1.0) # indirect doctest
     502            0.388621540699059
     503
     504        """
     505        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     506        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     507        prec = parent.prec() if hasattr(parent, 'prec') else 53
     508        if algorithm == 'mpmath':
     509            import mpmath
     510            from sage.libs.mpmath import utils as mpmath_utils
     511            return mpmath_utils.call(mpmath.airybi, x, derivative=alpha, parent=parent)
     512        elif algorithm == 'maxima':
     513            raise NotImplementedError("general case not available in maxima")
     514        else:
     515            raise ValueError("unknown algorithm")
     516
     517class FunctionAiryBiSimple(BuiltinFunction):
     518    def __init__(self):
     519        """
     520        The class for the Airy Bi function.
     521
     522        EXAMPLES::
     523
     524            sage: from sage.functions.airy import FunctionAiryBiSimple
     525            sage: airy_bi_simple= FunctionAiryBiSimple()
     526            sage: f=airy_bi_simple(x); f
     527            airy_bi(x)
     528        """
     529
     530        BuiltinFunction.__init__(self, "airy_bi", latex_name=r'\operatorname{Bi}',
     531                conversions=dict(mathematica='AiryBi', maxima='airy_bi'))
     532
     533    def _derivative_(self, x, diff_param=None):
     534        """
     535        EXAMPLES::
     536       
     537            sage: derivative(airy_bi(x),x) # indirect doctest
     538            airy_bi_prime(x)
     539        """
     540        return airy_bi_prime(x)
     541
     542    def _eval_(self, x):
     543        """
     544        EXAMPLES::
     545       
     546            sage: airy_bi(0) # indirect doctest
     547            1/3*3^(5/6)/gamma(1/3)
     548            sage: airy_bi(0.0) # indirect doctest
     549            0.614926627446001
     550            sage: airy_bi(I) # indirect doctest
     551            airy_bi(I)
     552            sage: airy_bi(1.0*I) # indirect doctest
     553            0.648858208330395 + 0.344958634768048*I
     554        """
     555
     556        if not isinstance(x,Expression):
     557            if is_inexact(x):
     558                return self._evalf_(x)
     559            elif x==0:
     560                one_sixth = ZZ(1)/6
     561                return 1/(3**(one_sixth)*gamma(2*one_sixth))
     562        else:
     563            return None
     564
     565    def _evalf_(self, x, **kwargs):
     566        """
     567        EXAMPLES::
     568       
     569            sage: airy_bi(0.0) # indirect doctest
     570            0.614926627446001
     571            sage: airy_bi(1.0*I) # indirect doctest
     572            0.648858208330395 + 0.344958634768048*I
     573
     574        We can use several methods for numerical evaluation::
     575
     576            sage: airy_bi(3).n(algorithm='maxima')
     577            14.0373289637
     578            sage: airy_bi(3).n(algorithm='mpmath')
     579            14.0373289637302
     580            sage: airy_bi(3).n(algorithm='mpmath', prec=100)
     581            14.037328963730232031740267314
     582
     583        TESTS::
     584
     585            sage: airy_bi(3).n(algorithm='maxima', prec=100)
     586            Traceback (most recent call last):
     587            ...
     588            ValueError: for the maxima algorithm the precision must be 53
     589
     590        """
     591        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     592        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     593        prec = parent.prec() if hasattr(parent, 'prec') else 53
     594        if algorithm == 'mpmath':
     595            import mpmath
     596            from sage.libs.mpmath import utils as mpmath_utils
     597            return mpmath_utils.call(mpmath.airybi, x, parent=parent)
     598        elif algorithm == 'maxima':
     599            if prec != 53:
     600                raise ValueError, "for the maxima algorithm the precision must be 53"
     601            return RDF(meval("airy_bi(%s)" % RDF(x)))
     602        else:
     603            raise ValueError("unknown algorithm")
     604
     605class FunctionAiryBiPrime(BuiltinFunction):
     606    def __init__(self):
     607        """
     608        The derivative of the Airy Bi function; see :func:`airy_bi`
     609        for the full documentation.
     610
     611        EXAMPLES::
     612
     613            sage: x,n=var('x n')
     614            sage: airy_bi_prime(x) # indirect doctest
     615            airy_bi_prime(x)
     616            sage: airy_bi_prime(0) # indirect doctest
     617            3^(1/6)/gamma(1/3)
     618        """
     619
     620        BuiltinFunction.__init__(self, "airy_bi_prime",
     621                latex_name=r"\operatorname{Bi}'",
     622                conversions=dict(mathematica='AiryBiPrime', maxima='airy_dbi'))
     623
     624    def _derivative_(self, x, diff_param=None):
     625        """
     626        EXAMPLES::
     627
     628            sage: derivative(airy_bi_prime(x),x) # indirect doctest
     629            x*airy_bi(x)
     630        """
     631        return x*airy_bi_simple(x)
     632
     633    def _eval_(self, x):
     634        """
     635        EXAMPLES::
     636
     637            sage: airy_bi(1,0) # indirect doctest
     638            3^(1/6)/gamma(1/3)
     639            sage: airy_bi(1,0.0) # indirect doctest
     640            0.448288357353826
     641        """
     642        if not isinstance(x,Expression):
     643            if is_inexact(x):
     644                return self._evalf_(x)
     645            elif x==0:
     646                one_sixth = ZZ(1)/6
     647                return 3**(one_sixth)/gamma(2*one_sixth)
     648        else:
     649            return None
     650
     651    def _evalf_(self, x, **kwargs):
     652        """
     653        EXAMPLES::
     654       
     655            sage: airy_bi(1,0.0) # indirect doctest
     656            0.448288357353826
     657
     658        We can use several methods for numerical evaluation::
     659
     660            sage: airy_bi(1,4).n(algorithm='maxima')
     661            161.926683505
     662            sage: airy_bi(1,4).n(algorithm='mpmath')
     663            161.926683504613
     664            sage: airy_bi(1,4).n(algorithm='mpmath', prec=100)
     665            161.92668350461340184309492429
     666
     667        TESTS::
     668
     669            sage: airy_bi(1,4).n(algorithm='maxima', prec=70)
     670            Traceback (most recent call last):
     671            ...
     672            ValueError: for the maxima algorithm the precision must be 53
     673
     674        """
     675        algorithm = kwargs.get('algorithm', None) or 'mpmath'
     676        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
     677        prec = parent.prec() if hasattr(parent, 'prec') else 53
     678        if algorithm == 'mpmath':
     679            import mpmath
     680            from sage.libs.mpmath import utils as mpmath_utils
     681            return mpmath_utils.call(mpmath.airybi, x, derivative=1, parent=parent)
     682        elif algorithm == 'maxima':
     683            if prec != 53:
     684                raise ValueError, "for the maxima algorithm the precision must be 53"
     685            return RDF(meval("airy_dbi(%s)" % RDF(x)))
     686        else:
     687            raise ValueError("unknown algorithm")
     688
     689airy_bi_general=FunctionAiryBiGeneral()
     690airy_bi_simple= FunctionAiryBiSimple()
     691airy_bi_prime=  FunctionAiryBiPrime()
     692
     693def airy_bi(alpha,x=None, hold_derivative=False, *args, **kwds):
     694    r"""
     695    The Airy Bi function `\operatorname{Bi}(x)` is one of the two
     696    linearly independent solutions to the Airy differental equation
     697    `f''(z) +f(z)x=0`, defined by the initial conditions:
     698
     699    .. math ::
     700        \operatorname{Bi}(0)=\frac{1}{3^{1/6} \Gamma(\frac{2}{3})},
     701
     702        \operatorname{Bi}'(0)=\frac{3^{1/6}}{ \Gamma(\frac{1}{3})}.
     703
     704    Another way to define the Airy Bi function is:
     705
     706    .. math::
     707        \operatorname{Bi}(x)=\frac{1}{\pi}\int_0^\infty
     708        \left[ \exp\left( xt -\frac{t^3}{3} \right)
     709        +\sin\left(xt + \frac{1}{3}t^3\right) \right ] dt.
     710
     711    INPUT:
     712   
     713    - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`.
     714      For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Bi}^{(n)}(z)`,
     715      and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
     716     
     717    .. math ::
     718     
     719        f_0(z) = \operatorname{Bi}(z)
     720     
     721        f_n(z) = \int_0^z f_{n-1}(t) dt.
     722
     723    - ``x``.- The argument of the function.
     724
     725    - ``hold_derivative``.- Whether or not to stop from returning higher derivatives in
     726      terms of `\operatorname{Bi}(x)` and `\operatorname{Bi}'(x)`.
     727   
     728
     729    EXAMPLES::
     730
     731        sage: n,x=var('n x')
     732        sage: airy_bi(x)
     733        airy_bi(x)
     734
     735    It can return derivatives or integrals::
     736   
     737        sage: airy_bi(1,x)
     738        airy_bi_prime(x)
     739        sage: airy_bi(2,x)
     740        x*airy_bi(x)
     741        sage: airy_bi(2,x,True)
     742        airy_bi(2, x)
     743        sage: airy_bi(-2,x)
     744        airy_bi(-2, x)
     745        sage: airy_bi(n, x)
     746        airy_bi(n, x)
     747
     748    It can be evaluated symbolically or numerically for real or complex values::
     749
     750        sage: airy_bi(0)
     751        1/3*3^(5/6)/gamma(1/3)
     752        sage: airy_bi(0.0)
     753        0.614926627446001
     754        sage: airy_bi(I)
     755        airy_bi(I)
     756        sage: airy_bi(1.0*I)
     757        0.648858208330395 + 0.344958634768048*I
     758
     759    The functions can be evaluated numerically using either mpmath (the default)
     760    or maxima, where mpmath can compute the values to arbitrary precision::
     761
     762        sage: airy_bi(2).n(prec=100)
     763        3.2980949999782147102806044252
     764        sage: airy_bi(2).n(algorithm='mpmath',prec=100)
     765        3.2980949999782147102806044252
     766        sage: airy_bi(2).n(algorithm='maxima')
     767        3.29809499998
     768
     769    And the derivatives can be evaluated::
     770
     771        sage: airy_bi(1,0)
     772        3^(1/6)/gamma(1/3)
     773        sage: airy_bi(1,0.0)
     774        0.448288357353826
     775
     776    Plots::
     777
     778        sage: plot(airy_bi(x),(x,-10,5))+plot(airy_bi_prime(x),(x,-10,5),color='red')
     779       
     780    **References**
     781     
     782    - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10"
     783
     784    - http://en.wikipedia.org/wiki/Airy_function
     785    """
     786
     787    #We catch the case with no alpha
     788    if x==None:
     789        x=alpha
     790        return airy_bi_simple(x, **kwds)
     791    #We raise an error if there are too many arguments
     792    if len(args) > 0:
     793        raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" %
     794                        (len(args)+3))
     795
     796    #We take care of all other cases.
     797    if hold_derivative:
     798        return airy_bi_general(alpha,x,**kwds)
     799    elif alpha==0:
     800        return airy_bi_simple(x, **kwds)
     801    elif alpha==1:
     802        return airy_bi_prime(x, **kwds)
     803    elif alpha>1:
     804        #we use a different variable here because if x is a
     805        #particular value, we would be differentiating a constant
     806        #which would return 0. What we want is the value of
     807        #the derivative at the value and not the derivative of
     808        #a particular value of the function.
     809        v=var('v')
     810        return derivative(airy_bi_simple(v,**kwds),v,alpha).subs(v=x)
     811    else:
     812        return airy_bi_general(alpha,x,**kwds)
     813
  • sage/functions/all.py

    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    3535                     inverse_jacobi,
    3636                     lngamma, error_fcn, elliptic_e,
    3737                     elliptic_f, elliptic_ec, elliptic_eu,
    38                      elliptic_kc, elliptic_pi, elliptic_j,
    39                      airy_ai, airy_bi)
     38                     elliptic_kc, elliptic_pi, elliptic_j)
    4039
    4140from orthogonal_polys import (chebyshev_T,
    4241                              chebyshev_U,
     
    6362
    6463from min_max import max_symbolic, min_symbolic
    6564
     65from airy import airy_ai, airy_ai_prime, airy_bi, airy_bi_prime
     66
    6667from exp_integral import (exp_integral_e, exp_integral_e1, log_integral, li, Li,
    6768                          log_integral_offset,
    6869                          sin_integral, cos_integral, Si, Ci,
  • sage/functions/special.py

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    3030implemented here.
    3131
    3232
    33 -  Airy function The function `Ai(x)` and the related
    34    function `Bi(x)`, which is also called an Airy function,
    35    are solutions to the differential equation
    36 
    37    
    38    .. math::
    39 
    40          y'' - xy = 0,
    41 
    42    known as the Airy equation. They belong to the class of 'Bessel functions of
    43    fractional order'. The initial conditions
    44    `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
    45    `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define
    46    `Ai(x)`. The initial conditions
    47    `Bi(0) = 3^{1/2}Ai(0)`, `Bi'(0) = -3^{1/2}Ai'(0)`
    48    define `Bi(x)`.
    49 
    50    They are named after the British astronomer George Biddell Airy.
    51 
    5233-  Spherical harmonics: Laplace's equation in spherical coordinates
    5334   is:
    5435   
     
    291272- Abramowitz and Stegun: Handbook of Mathematical Functions,
    292273  http://www.math.sfu.ca/~cbm/aands/
    293274
    294 - http://en.wikipedia.org/wiki/Airy_function
    295 
    296275- http://en.wikipedia.org/wiki/Spherical_harmonics
    297276
    298277- http://en.wikipedia.org/wiki/Helmholtz_equation
     
    362341
    363342    Then after using one of these functions, it changes::
    364343
    365         sage: from sage.functions.special import airy_ai
    366         sage: airy_ai(1.0)
    367         0.135292416313
     344        sage: from sage.functions.special import spherical_harmonic
     345        sage: spherical_harmonic(3,2,1,2)
     346        15/4*sqrt(7/30)*cos(1)*e^(4*I)*sin(1)^2/sqrt(pi)
    368347        sage: sage.functions.special._done
    369348        True
    370349    """
     
    382361   
    383362    TEST::
    384363
    385         sage: from sage.functions.special import airy_ai
    386         sage: airy_bi(1.0)
    387         1.20742359495
     364        sage: from sage.functions.special import spherical_bessel_J
     365        sage: spherical_bessel_J(2.,3.)
     366        0.298637497076
     367
    388368    """
    389369    return maxima(x).sage()
    390370
     
    544524    return NewMaximaFunction()
    545525
    546526
    547 def airy_ai(x):
    548    r"""
    549    The function `Ai(x)` and the related function `Bi(x)`,
    550    which is also called an *Airy function*, are
    551    solutions to the differential equation
    552 
    553    .. math::
    554 
    555       y'' - xy = 0,
    556 
    557    known as the *Airy equation*. The initial conditions
    558    `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
    559    `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define `Ai(x)`.
    560    The initial conditions `Bi(0) = 3^{1/2}Ai(0)`,
    561    `Bi'(0) = -3^{1/2}Ai'(0)` define `Bi(x)`.
    562 
    563    They are named after the British astronomer George Biddell Airy.
    564    They belong to the class of "Bessel functions of fractional order".
    565 
    566    EXAMPLES::
    567    
    568        sage: airy_ai(1.0)        # last few digits are random
    569        0.135292416312881400
    570        sage: airy_bi(1.0)        # last few digits are random
    571        1.20742359495287099
    572 
    573    REFERENCE:
    574 
    575    - Abramowitz and Stegun: Handbook of Mathematical Functions,
    576      http://www.math.sfu.ca/~cbm/aands/
    577 
    578    - http://en.wikipedia.org/wiki/Airy_function
    579    """
    580    _init()
    581    return RDF(meval("airy_ai(%s)"%RDF(x)))
    582 
    583 def airy_bi(x):
    584    r"""
    585    The function `Ai(x)` and the related function `Bi(x)`,
    586    which is also called an *Airy function*, are
    587    solutions to the differential equation
    588 
    589    .. math::
    590 
    591       y'' - xy = 0,
    592 
    593    known as the *Airy equation*. The initial conditions
    594    `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
    595    `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define `Ai(x)`.
    596    The initial conditions `Bi(0) = 3^{1/2}Ai(0)`,
    597    `Bi'(0) = -3^{1/2}Ai'(0)` define `Bi(x)`.
    598 
    599    They are named after the British astronomer George Biddell Airy.
    600    They belong to the class of "Bessel functions of fractional order".
    601 
    602    EXAMPLES::
    603    
    604        sage: airy_ai(1)        # last few digits are random
    605        0.135292416312881400
    606        sage: airy_bi(1)        # last few digits are random
    607        1.20742359495287099
    608 
    609    REFERENCE:
    610 
    611    - Abramowitz and Stegun: Handbook of Mathematical Functions,
    612      http://www.math.sfu.ca/~cbm/aands/
    613 
    614    - http://en.wikipedia.org/wiki/Airy_function
    615    """
    616    _init()
    617    return RDF(meval("airy_bi(%s)"%RDF(x)))
    618 
    619 
    620527def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53):
    621528    r"""
    622529    Default is a wrap of PARI's hyperu(alpha,beta,x) function.