Search code examples
pythonscipydata-analysisphysics

E1/3 function in python?

I am analyzing some lab results with python and I need to fit the data to a distribution that includes the generalized exponential integral E1/3. The closest thing I could find besides calculating the integral myself (which is very slow) is scipy.special.expn(n, x), but it only supports integer values for n, and I need n=1/3.

Does anybody know if such a function exists?

Solution

  • Updated

    If writing your own using a quadrature routine isn't accurate enough you can consider using sympy. The function is implemented in sympy.functions.special.error_functions.expint. Since you want a numerical solution, you'll need to convert it to a float.

    import sympy as sp

def expn(n, x):
    return float(sp.functions.special.error_functions.expint(n, x))


print(expn(1/3,1))  # 0.3044294477841587

    It's not fast, but it might do what you want.

    Edit: Another alternative is to use the definition of the exponential integral based on the upper incomplete gamma function (scipy has a normalized function for this), as per the digital library of mathematical functions.

    from scipy.special import gammaincc, gamma

def expn(n, x):
    return x**(n-1)*gammaincc(1-n,x)*gamma(1-n)

print(expn(1/3, 1)) # 0.3044294477841586

    Old Solution

    Osama's solution using quadrature isn't general enough, so I am including my own version below.

    from scipy.integrate import quad
import numpy as np

def expn(n, x):
    return quad(lambda t: np.exp(-x*t)/t**n, 1, np.inf)[0]

print(expn(1/3,1))  # 0.3044294477840883

    In terms of accuracy, Mathematica reports the result to be 0.30442944778415870046. So, the above quadrature solution has a 2.312e-11% error. I can't imagine a higher accuracy solution would be necessary.