Ticket #11143: trac_11143_En.patch

File trac_11143_En.patch, 7.3 KB (added by benjaminfjones, 10 years ago)

added the symbolic function exp_integral_e

  • sage/functions/all.py

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1308594346 14400
    # Node ID 777dfa0aacecef5e9db8d6552c9e832c96d599bf
    # Parent  20a5478d27a7ba3a6662977094212ece3d187a7b
    Trac 11143: added symbolic function `exp_integral_e`
    
    diff -r 20a5478d27a7 -r 777dfa0aacec sage/functions/all.py
    a b  
    3838                     lngamma, exp_int, error_fcn, elliptic_e,
    3939                     elliptic_f, elliptic_ec, elliptic_eu,
    4040                     elliptic_kc, elliptic_pi, elliptic_j,
    41                      airy_ai, airy_bi)
     41                     airy_ai, airy_bi,
     42                     exp_integral_e)
    4243                       
    4344from orthogonal_polys import (chebyshev_T,
    4445                              chebyshev_U,
  • sage/functions/special.py

    diff -r 20a5478d27a7 -r 777dfa0aacec sage/functions/special.py
    a b  
    1212
    1313- David Joyner (2008-04-23): addition of elliptic integrals
    1414
     15- Benjamin Jones (2011-06-12): addition of the generalized exponential integral
     16
    1517This module provides easy access to many of Maxima and PARI's
    1618special functions.
    1719
     
    393395import sage.rings.commutative_ring as commutative_ring
    394396import sage.rings.ring as ring
    395397from sage.functions.other import real, imag
    396 from sage.symbolic.function import BuiltinFunction
     398from sage.functions.log import exp
     399from sage.symbolic.function import BuiltinFunction, is_inexact
    397400from sage.calculus.calculus import maxima
     401from sage.symbolic.expression import Expression
     402from sage.structure.parent import Parent
     403from sage.structure.coerce import parent
     404from sage.libs.mpmath import utils as mpmath_utils
    398405
    399406_done = False
    400407def _init():
     
    17741781        except:
    17751782            raise NotImplementedError
    17761783
     1784class Function_exp_integral_e(BuiltinFunction):
     1785    r"""
     1786    The generalized complex exponential integral `E_n(z)` defined by
    17771787
     1788    .. math::
    17781789
     1790        \operatorname{E_n}(z) = \int_1^{\infty} \frac{e^{-z t}}{t^n} \; dt
     1791
     1792    for complex numbers `n` and `z`, see [AS]_ 5.1.4.
     1793
     1794    The special case where `n = 1` is denoted in Sage by
     1795    :meth:`exponential_integral_1`.
     1796
     1797    EXAMPLES:
     1798
     1799    Numerical evaluation is handled using mpmath::
     1800
     1801        sage: N(exp_integral_e(1,1))
     1802        0.219383934395520
     1803        sage: exp_integral_e(1, RealField(100)(1))
     1804        0.21938393439552027367716377546
     1805
     1806    We can compare this to PARI's evaluation of
     1807    :meth:`exponential_integral_1`::
     1808
     1809        sage: N(exponential_integral_1(1))
     1810        0.219383934395520
     1811
     1812    We can verify one case of [AS]_ 5.1.45, i.e.
     1813    `E_n(z) = z^{n-1}\Gamma(1-n,z)`::
     1814
     1815        sage: N(exp_integral_e(2, 3+I))
     1816        0.00354575823814662 - 0.00973200528288687*I
     1817        sage: N((3+I)*gamma(-1, 3+I))
     1818        0.00354575823814662 - 0.00973200528288687*I
     1819
     1820    Maxima returns the following improper integral as a multiple of
     1821    exp_integral_e(1,1)::
     1822
     1823        sage: uu = integral(e^(-x)*log(x+1),x,0,oo)
     1824        sage: uu
     1825        e*exp_integral_e(1, 1)
     1826        sage: uu.n(digits=30)
     1827        0.596347362323194074341078499369
     1828
     1829    Symbolic derivatives and integrals are handled by Sage and Maxima::
     1830
     1831        sage: x = var('x')
     1832        sage: f = exp_integral_e(2,x)
     1833        sage: f.diff(x)
     1834        -exp_integral_e(1, x)
     1835       
     1836        sage: f.integrate(x)
     1837        -exp_integral_e(3, x)
     1838       
     1839        sage: f = exp_integral_e(-1,x)
     1840        sage: f.integrate(x)
     1841        Ei(-x) - gamma(-1, x)
     1842       
     1843    Some special values of `exp_integral_e` can be simplified.
     1844    [AS]_ 5.1.23::
     1845       
     1846        sage: f = exp_integral_e(0,x)
     1847        sage: f.simplify()
     1848        e^(-x)/x
     1849       
     1850    [AS]_ 5.1.24::
     1851   
     1852        sage: nn = var('nn')
     1853        sage: assume(nn > 1)
     1854        sage: f = exp_integral_e(nn,0)
     1855        sage: f.simplify()
     1856        1/(nn - 1)
     1857
     1858    ALGORITHM:
     1859
     1860    Numerical evaluation is handled using mpmath, but symbolics are handled
     1861    by Sage and Maxima.
     1862
     1863    REFERENCES:
     1864
     1865    .. _[AS] 'Handbook of Mathematical Functions', Milton Abramowitz and Irene
     1866    A. Stegun, National Bureau of Standards Applied Mathematics Series, 55.
     1867    See also http://www.math.sfu.ca/~cbm/aands/.
     1868
     1869    """
     1870    def __init__(self):
     1871        """
     1872        See the docstring for :class:`Function_exp_integral_e`.
     1873
     1874        EXAMPLES::
     1875
     1876            sage: exp_integral_e(1,0)
     1877            exp_integral_e(1, 0)
     1878
     1879        """
     1880        BuiltinFunction.__init__(self, "exp_integral_e", nargs=2, latex_name=r'exp_integral_e',
     1881                                 conversions=dict(maxima='expintegral_e'))
     1882
     1883    def _eval_(self, n, z):
     1884        """
     1885        EXAMPLES::
     1886
     1887            sage: exp_integral_e(1.0, x)
     1888            exp_integral_e(1.00000000000000, x)
     1889            sage: exp_integral_e(x, 1.0)
     1890            exp_integral_e(x, 1.00000000000000)
     1891            sage: exp_integral_e(1.0, 1.0)
     1892            0.219383934395520
     1893           
     1894        """
     1895        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
     1896               (is_inexact(n) or is_inexact(z)):
     1897            coercion_model = sage.structure.element.get_coercion_model()
     1898            n, z = coercion_model.canonical_coercion(n, z)
     1899            return self._evalf_(n, z, parent(n))
     1900       
     1901        z_zero = False   
     1902        # special case: z == 0 and n > 1
     1903        if isinstance(z, Expression):
     1904            if z._is_numerically_zero():
     1905                z_zero = True # for later
     1906                if n > 1:
     1907                    return 1/(n-1)
     1908        else:
     1909            if not z:
     1910                z_zero = True
     1911                if n > 1:
     1912                    return 1/(n-1)
     1913
     1914        # special case: n == 0
     1915        if isinstance(n, Expression):
     1916            if n._is_numerically_zero():
     1917                if z_zero:
     1918                    return None
     1919                else:
     1920                    return exp(-z)/z
     1921        else:
     1922            if not n:
     1923                if z_zero:
     1924                    return None
     1925                else:
     1926                    return exp(-z)/z
     1927
     1928        return None # leaves the expression unevaluated
     1929
     1930    def _evalf_(self, n, z, parent=None):
     1931        """
     1932        EXAMPLES::
     1933
     1934            sage: N(exp_integral_e(1, 1+I))
     1935            0.000281624451981418 - 0.179324535039359*I
     1936            sage: exp_integral_e(1, RealField(100)(1))
     1937            0.21938393439552027367716377546
     1938
     1939        """
     1940        import mpmath
     1941        return mpmath_utils.call(mpmath.expint, n, z, parent=parent)
     1942       
     1943    def _derivative_(self, n, z, diff_param=None):
     1944        """
     1945        If `n` is an integer strictly larger than 0, then the derivative of
     1946        `exp_integral_e(n,z)` with respect to `z` is
     1947        `-1*exp_integral_e(n-1,z)`. See [AS], 5.1.26.
     1948
     1949        EXAMPLES::
     1950
     1951            sage: x = var('x')
     1952            sage: f = exp_integral_e(2,x)
     1953            sage: f.diff(x)
     1954            -exp_integral_e(1, x)
     1955
     1956            sage: f = exp_integral_e(2,sqrt(x))
     1957            sage: f.diff(x)
     1958            -1/2*exp_integral_e(1, sqrt(x))/sqrt(x)
     1959
     1960        """
     1961        if n in ZZ and n > 0:
     1962            return -1*exp_integral_e(n-1,z)
     1963        else:
     1964            raise NotImplementedError("The derivative of is only implemented for n = 1, 2, 3, ...")
     1965
     1966exp_integral_e = Function_exp_integral_e()