更新

如果使用正交例程编写自己的程序不够准确，您可以考虑使用sympy。该函数在sympy.functions.special.error_functions.expint中实现。因为你想要一个数值解，你需要把它转换成一个浮点数。

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

它不快，但它可能会做你想要的。
编辑：另一种选择是使用基于上不完全伽马函数的指数积分的定义（scipy对此有一个归一化函数），根据数学函数的数字库。

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

旧方案

奥萨马使用正交的解决方案不够通用，所以我在下面包括我自己的版本。

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

在精度方面，Mathematica报告结果为0.30442944778415870046。因此，上述正交解具有2.312e-11%的误差。我无法想象一个更高精度的解决方案将是必要的。

我认为你需要定义函数来计算这个，

import numpy as np
from scipy.integrate import quad

def integrand(x):
    return np.exp(-x) / x**(1/3)

result, error = quad(integrand, 0, np.inf)
print(result)

但它适用于学术用途，我的导师推荐了一个python库“mpmath”ref：offical mpmath website，比上面的函数快得多。所以修改后的snippet将是这样的：

import mpmath

def integrand(x):
    return mpmath.exp(-x) / x**(1/3)

result = mpmath.quad(integrand, [0, mpmath.inf])
print(result)