Ticket #14896: trac_14896_2.patch

File trac_14896_2.patch, 13.9 KB (added by eviatarbach, 9 years ago)

  • sage/functions/all.py 

    # HG changeset patch
# User Eviatar Bach <eviatarbach@gmail.com>
# Date 1374019480 25200
# Branch hyper
# Node ID 93aa2f7392e1909a85a3571f3d93dc54f0996083
# Parent  2c773352ebc9bc079a27fa06e03ca858bf60e39a
Implementing confluent hypergeometric functions

diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    2727from transcendental import (zeta, zetaderiv, zeta_symmetric, dickman_rho)
    2828
    2929from special import (bessel_I, bessel_J, bessel_K, bessel_Y,
    30                      hypergeometric_U, Bessel,
     30                     Bessel,
    3131                     spherical_bessel_J, spherical_bessel_Y,
    3232                     spherical_hankel1, spherical_hankel2,
    3333                     spherical_harmonic, jacobi,
    … …  
    6868                          sinh_integral, cosh_integral, Shi, Chi,
    6969                          exponential_integral_1, Ei)
    7070
    71 from hypergeometric import hypergeometric
     71from hypergeometric import hypergeometric, hypergeometric_M, hypergeometric_U

  • sage/functions/hypergeometric.py 

    diff --git a/sage/functions/hypergeometric.py b/sage/functions/hypergeometric.py
    a b  
    11r"""
    22Hypergeometric Functions
    33
    4 This module implements manipulation of infinite hypergeometric series
    5 represented in standard parametric form (as $\,_pF_q$ functions).
     4This module implements manipulation of generalized hypergeometric series
     5represented in standard parametric form (as $\,_pF_q$ functions), as well as
     6the confluent hypergeometric functions of the first and second kind.
    67
    78AUTHORS:
    89
    … …  
    123124    x),), (3,), x)
    124125    sage: nest(lambda y: hypergeometric([y], [], x), 3, 1)._mathematica_init_()
    125126    'HypergeometricPFQ[{HypergeometricPFQ[{HypergeometricPFQ[{1},{},x]},...
     127
     128The confluent hypergeometric functions can arise as solutions to second-order
     129differential equations (example from `here <http://ask.sagemath.org/question/
     1301168/how-can-one-use-maxima-kummer-confluent-functions>`_)::
     131
     132    sage: m = var('m')
     133    sage: y = function('y', x)
     134    sage: desolve(diff(y, x, 2) + 2*x*diff(y, x) - 4*m*y, y,
     135    ....:         contrib_ode=true, ivar=x)
     136    [y(x) == k1*hypergeometric_M(-m, 1/2, -x^2) +...
     137     k2*hypergeometric_U(-m, 1/2, -x^2)]
     138
     139Series expansions of confluent hypergeometric functions::
     140
     141    sage: hypergeometric_M(2, 2, x).series(x, 3)
     142    1 + 1*x + 1/2*x^2 + Order(x^3)
     143    sage: hypergeometric_U(2, 2, x).series(x == 3, 100).subs(x=1).n()
     144    0.403652637676806
     145    sage: hypergeometric_U(2, 2, 1).n()                             
     146    0.403652637676806
    126147"""
    127148#*****************************************************************************
    128149#       Copyright (C) 2010 Fredrik Johansson <fredrik.johansson@gmail.com>
    … …  
    138159                            rising_factorial, factorial)
    139160from sage.functions.other import sqrt, gamma, real_part
    140161from sage.functions.log import exp, log
    141 from sage.functions.trig import cos, sin
     162from sage.functions.trig import sin
    142163from sage.functions.hyperbolic import cosh, sinh
    143164from sage.functions.other import erf
    144165from sage.symbolic.constants import pi
    … …  
    782803
    783804            def _2f1(a, b, c, z):
    784805                """
    785                 Evaluation of 2F1(a, b, c, z), assuming a, b, c positive
     806                Evaluation of 2F1(a, b; c; z), assuming a, b, c positive
    786807                integers or half-integers
    787808                """
    788809                if b == c:
    … …  
    838859                    pass
    839860        return hyp
    840861    return sum([coeff * _closed_form(pfq) for coeff, pfq in new._deflated()])
     862
     863class Hypergeometric_M(BuiltinFunction):
     864    r"""
     865    The confluent hypergeometric function of the first kind,
     866    `y = M(a,b,z)`, is defined to be the solution to Kummer's differential
     867    equation
     868   
     869    .. math::
     870   
     871             zy'' + (b-z)y' - ay = 0.
     872
     873    This is not the same as Kummer's `U`-hypergeometric function, though it
     874    satisfies the same DE that `M` does.
     875   
     876    .. warning::
     877
     878       In the literature, both are called "Kummer confluent
     879       hypergeometric" functions.
     880   
     881    EXAMPLES::
     882   
     883        sage: hypergeometric_M(1, 1, 1)
     884        hypergeometric_M(1, 1, 1)
     885        sage: hypergeometric_M(1, 1, 1.)
     886        2.71828182845905
     887        sage: hypergeometric_M(1, 1, 1).n(70)
     888        2.7182818284590452354
     889        sage: hypergeometric_M(1, 1, 1).simplify_hypergeometric()
     890        e
     891        sage: hypergeometric_M(1, 1/2, x).simplify_hypergeometric()
     892        (-I*sqrt(pi)*x*e^x*erf(I*sqrt(-x)) + sqrt(-x))/sqrt(-x)
     893        sage: hypergeometric_M(1, 3/2, 1).simplify_hypergeometric()
     894        1/2*sqrt(pi)*e*erf(1)
     895    """
     896    def __init__(self):
     897        BuiltinFunction.__init__(self, 'hypergeometric_M', nargs=3,
     898                                 conversions={'mathematica':
     899                                              'Hypergeometric1F1',
     900                                              'maxima': 'kummer_m'},
     901                                 latex_name='M')
     902
     903    def _eval_(self, a, b, z, **kwargs):
     904        cm = get_coercion_model()
     905        co = cm.canonical_coercion(cm.canonical_coercion(a, b)[0], z)[0]
     906        if is_inexact(co) and not isinstance(co, Expression):
     907            from sage.structure.coerce import parent
     908            return self._evalf_(a, b, z, parent=parent(co))
     909        if not isinstance(z, Expression) and z == 0:
     910            return Integer(1)
     911        return
     912
     913    def _evalf_(self, a, b, z, parent):
     914        from mpmath import hyp1f1
     915        return mpmath_utils.call(hyp1f1, a, b, z, parent=parent)
     916
     917    def _derivative_(self, a, b, z, diff_param):
     918        if diff_param == 2:
     919            return (a / b) * hypergeometric_M(a + 1, b + 1, z)
     920        raise NotImplementedError('derivative of hypergeometric function '
     921                                  'with respect to parameters')
     922
     923    class EvaluationMethods:
     924        def generalized(cls, self, a, b, z):
     925            """
     926            Return as a generalized hypergeometric function
     927
     928            EXAMPLES::
     929
     930                sage: a, b, z = var('a b z')
     931                sage: hypergeometric_M(a, b, z).generalized()
     932                hypergeometric((a,), (b,), z)
     933
     934            """
     935            return hypergeometric([a], [b], z)
     936
     937hypergeometric_M = Hypergeometric_M()
     938
     939class Hypergeometric_U(BuiltinFunction):
     940    r"""
     941    The confluent hypergeometric function of the second kind,
     942    `y = U(a,b,z)`, is defined to be the solution to Kummer's differential
     943    equation
     944   
     945    .. math::
     946   
     947             zy'' + (b-z)y' - ay = 0.
     948
     949    This satisfies `U(a,b,z) \sim z^{-a}`, as
     950    `z\rightarrow \infty`, and is sometimes denoted
     951    `z^{-a}{}_2F_0(a,1+a-b;;-1/z)`. This is not the same as Kummer's
     952    `M`-hypergeometric function, denoted sometimes as
     953    `_1F_1(\alpha,\beta,z)`, though it satisfies the same DE that
     954    `U` does.
     955   
     956    .. warning::
     957
     958       In the literature, both are called "Kummer confluent
     959       hypergeometric" functions.
     960   
     961    EXAMPLES::
     962   
     963        sage: hypergeometric_U(1, 1, 1)
     964        hypergeometric_U(1, 1, 1)
     965        sage: hypergeometric_U(1, 1, 1.)
     966        0.596347362323194
     967        sage: hypergeometric_U(1, 1, 1).n(70)
     968        0.59634736232319407434
     969        sage: hypergeometric_U(10^4, 1/3, 1).n()
     970        6.60377008885811e-35745
     971        sage: hypergeometric_U(2 + I, 2, 1).n()
     972        0.183481989942099 - 0.458685959185190*I
     973        sage: hypergeometric_U(1, 3, x).simplify_hypergeometric()
     974        (x + 1)/x^2
     975        sage: hypergeometric_U(1, 2, 2).simplify_hypergeometric()
     976        1/2
     977
     978    """
     979    def __init__(self):
     980        BuiltinFunction.__init__(self, 'hypergeometric_U', nargs=3,
     981                                 conversions={'mathematica':
     982                                              'HypergeometricU',
     983                                              'maxima': 'kummer_u'},
     984                                 latex_name='U')
     985
     986    def _eval_(self, a, b, z, **kwargs):
     987        cm = get_coercion_model()
     988        co = cm.canonical_coercion(cm.canonical_coercion(a, b)[0], z)[0]
     989        if is_inexact(co) and not isinstance(co, Expression):
     990            from sage.structure.coerce import parent
     991            return self._evalf_(a, b, z, parent=parent(co))
     992        return
     993
     994    def _evalf_(self, a, b, z, parent):
     995        from mpmath import hyperu
     996        return mpmath_utils.call(hyperu, a, b, z, parent=parent)
     997
     998    def _derivative_(self, a, b, z, diff_param):
     999        if diff_param == 2:
     1000            return -a * hypergeometric_U(a + 1, b + 1, z)
     1001        raise NotImplementedError('derivative of hypergeometric function '
     1002                                  'with respect to parameters')
     1003
     1004    class EvaluationMethods:
     1005        def generalized(cls, self, a, b, z):
     1006            """
     1007            Return in terms of the generalized hypergeometric function
     1008
     1009            EXAMPLES::
     1010
     1011                sage: a, b, z = var('a b z')
     1012                sage: hypergeometric_U(a, b, z).generalized()
     1013                z^(-a)*hypergeometric((a, a - b + 1), (), -1/z)
     1014                sage: hypergeometric_U(1, 3, 1/2).generalized()           
     1015                2*hypergeometric((1, -1), (), -2)
     1016                sage: hypergeometric_U(3, I, 2).generalized() 
     1017                1/8*hypergeometric((3, -I + 4), (), -1/2)
     1018
     1019            """
     1020            return z ** (-a) * hypergeometric([a, a - b + 1], [], -z ** (-1))
     1021
     1022hypergeometric_U = Hypergeometric_U()

  • sage/functions/special.py 

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    175175         y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x)     = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x). 
    176176
    177177
    178 
    179 -  For `x>0`, the confluent hypergeometric function
    180    `y = U(a,b,x)` is defined to be the solution to Kummer's
    181    differential equation
    182 
    183    
    184    .. math::
    185 
    186      xy'' + (b-x)y' - ay = 0,
    187 
    188    which satisfies `U(a,b,x) \sim x^{-a}`, as
    189    `x\rightarrow \infty`. (There is a linearly independent
    190    solution, called Kummer's function `M(a,b,x)`, which is not
    191    implemented.)
    192 
    193178-  Jacobi elliptic functions can be thought of as generalizations
    194179   of both ordinary and hyperbolic trig functions. There are twelve
    195180   Jacobian elliptic functions. Each of the twelve corresponds to an
    … …  
    349334- Online Encyclopedia of Special Function
    350335  http://algo.inria.fr/esf/index.html
    351336
    352 TODO: Resolve weird bug in commented out code in hypergeometric_U
    353 below.
    354 
    355337AUTHORS:
    356338
    357339- David Joyner and William Stein
    … …  
    12161198            f = lambda z: bessel_Y(nu,z,s)
    12171199            P = plot(f,a,b)
    12181200        return P
    1219    
    1220 def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53):
    1221     r"""
    1222     Default is a wrap of PARI's hyperu(alpha,beta,x) function.
    1223     Optionally, algorithm = "scipy" can be used.
    1224    
    1225     The confluent hypergeometric function `y = U(a,b,x)` is
    1226     defined to be the solution to Kummer's differential equation
    1227    
    1228     .. math::
    1229    
    1230              xy'' + (b-x)y' - ay = 0.     
    1231    
    1232     This satisfies `U(a,b,x) \sim x^{-a}`, as
    1233     `x\rightarrow \infty`, and is sometimes denoted
    1234     ``x^{-a}2_F_0(a,1+a-b,-1/x)``. This is not the same as Kummer's
    1235     `M`-hypergeometric function, denoted sometimes as
    1236     ``_1F_1(alpha,beta,x)``, though it satisfies the same DE that
    1237     `U` does.
    1238    
    1239     .. warning::
    1240 
    1241        In the literature, both are called "Kummer confluent
    1242        hypergeometric" functions.
    1243    
    1244     EXAMPLES::
    1245    
    1246         sage: hypergeometric_U(1,1,1,"scipy")
    1247         0.596347362323...
    1248         sage: hypergeometric_U(1,1,1) 
    1249         0.59634736232319...
    1250         sage: hypergeometric_U(1,1,1,"pari",70)
    1251         0.59634736232319407434...
    1252     """
    1253     if algorithm=="scipy":
    1254         if prec != 53:
    1255             raise ValueError, "for the scipy algorithm the precision must be 53"
    1256         import scipy.special
    1257         ans = str(scipy.special.hyperu(float(alpha),float(beta),float(x)))
    1258         ans = ans.replace("(","")
    1259         ans = ans.replace(")","")
    1260         ans = ans.replace("j","*I")
    1261         return sage_eval(ans)
    1262     elif algorithm=='pari':
    1263         from sage.libs.pari.all import pari
    1264         R = RealField(prec)
    1265         return R(pari(R(alpha)).hyperu(R(beta), R(x), precision=prec))
    1266     else:
    1267         raise ValueError, "unknown algorithm '%s'"%algorithm
    12681201
    12691202def spherical_bessel_J(n, var, algorithm="maxima"):
    12701203    r"""

  • sage/symbolic/expression.pyx 

    diff --git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
    a b  
    74817481
    74827482    def simplify_hypergeometric(self, algorithm='maxima'):
    74837483        """
    7484         Simplify an expression containing hypergeometric functions
     7484        Simplify an expression containing hypergeometric or confluent
     7485        hypergeometric functions
    74857486
    74867487        INPUT:
    74877488
    … …  
    75097510            sage: (nest(lambda y: hypergeometric([y], [1], x), 3, 1)
    75107511            ....:  .simplify_hypergeometric(algorithm='sage'))
    75117512            hypergeometric((hypergeometric((e^x,), (1,), x),), (1,), x)
    7512 
    7513         """
    7514         from sage.functions.hypergeometric import hypergeometric, closed_form
     7513            sage: hypergeometric_M(1, 3, x).simplify_hypergeometric()     
     7514            -2*(x - e^x + 1)/x^2
     7515            sage: (2 * hypergeometric_U(1, 3, x)).simplify_hypergeometric()
     7516            2*(x + 1)/x^2
     7517
     7518        """
     7519        from sage.functions.hypergeometric import (hypergeometric,
     7520                                                   hypergeometric_M,
     7521                                                   hypergeometric_U,
     7522                                                   closed_form)
    75157523        from sage.calculus.calculus import maxima
    75167524        try:
    75177525            op = self.operator()
    75187526        except RuntimeError:
    75197527            return self
    75207528        ops = self.operands()
     7529        if op == hypergeometric_M or op == hypergeometric_U:
     7530            return self.generalized().simplify_hypergeometric(algorithm)
    75217531        if op == hypergeometric:
    75227532            if algorithm == 'maxima':
    75237533                return (self.parent()