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.