Python implementation of the Exponential Integral function¶
This small notebook is a Python 3 implementation of the Exponential Integral function, $Ei(x)$, defined like this:
$$ \forall x\in\mathbb{R}\setminus\{0\},\;\; \mathrm{Ei}(x) := \int_{-\infty}^x \frac{\mathrm{e}^u}{u} \; \mathrm{d}u. $$
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context="notebook", style="darkgrid", palette="hls", font="sans-serif", font_scale=1.4)
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (19.80, 10.80)
Two implementations¶
As one can show mathematically, there is another equivalent definition (the one used on Wikipedia):
$$ \forall x\in\mathbb{R}\setminus\{0\},\;\; \mathrm{Ei}(x) := - \int_{-x}^{\infty} \frac{\mathrm{e}^{-t}}{t} \; \mathrm{d}t. $$
Numerically, we will avoid the issue in $0$ by integrating up-to $-\varepsilon$ instead of $0^-$ and from $\varepsilon$ instead of $0^+$, for a some small $\varepsilon$ (e.g., $\varepsilon=10^{-7}$), and from $-M$ for a large value $M\in\mathbb{R}$ (e.g., $M=10000$), instead of $-\infty$.
We use the scipy.integrate.quad function.
from scipy.integrate import quad # need only 1 function
First definition is the simplest:
@np.vectorize
def Ei(x, minfloat=1e-7, maxfloat=10000):
"""Ei integral function."""
minfloat = min(np.abs(x), minfloat)
maxfloat = max(np.abs(x), maxfloat)
def f(t):
return np.exp(t) / t
if x > 0:
return (quad(f, -maxfloat, -minfloat)[0] + quad(f, minfloat, x)[0])
else:
return quad(f, -maxfloat, x)[0]
The other definition is very similar:
@np.vectorize
def Ei_2(x, minfloat=1e-7, maxfloat=10000):
"""Ei integral function."""
minfloat = min(np.abs(x), minfloat)
maxfloat = max(np.abs(x), maxfloat)
def f(t):
return np.exp(-t) / t
if x > 0:
return -1.0 * (quad(f, -x, -minfloat)[0] + quad(f, minfloat, maxfloat)[0])
else:
return -1.0 * quad(f, -x, maxfloat)[0]
Checking the two versions¶
We can quickly check that the two are equal:
from numpy.linalg import norm
X = np.linspace(-1, 1, 1000) # 1000 points
Y = Ei(X)
Y_2 = Ei_2(X)
assert np.allclose(Y, Y_2)
print(f"Two versions of Ei(x) are indeed equal for {len(X)} values.")
We can compare which is fastest to evaluate:
%timeit Y = Ei(X)
%timeit Y_2 = Ei_2(X)
They both take about the same time, but the second implementation seems (slightly) faster.
Comparison with scipy.special.expi¶
The $\mathrm{Ei}$ function is also implemented as scipy.special.expi:
from scipy.special import expi
Y_3 = expi(X)
np.allclose(Y, Y_3)
The difference is not too large:
np.max(np.abs(Y - Y_3))
assert np.allclose(Y, Y_3, rtol=1e-6, atol=1e-6)
print(f"Our version of Ei(x) is the same as the one in scipy.special.expi ({len(X)} values).")
Special values¶
We can compute some special values, like $\mathrm{Ei}(1)$ and solving (numerically) $\mathrm{Ei}(x)=0$.
Ei(1)
from scipy.optimize import root
res = root(Ei, x0=1)
res
print(f"The approximate solution to Ei(x)=0 is x0 = {res.x[0]} (for which Ei(x)={res.fun})...")
Limits¶
We can check that $\mathrm{Ei}(x)\to0$ for $x\to-\infty$ and $\mathrm{Ei}(x)\to+\infty$ for $x\to\infty$:
for x in -np.linspace(1, 1000, 10):
print(f"For x = {x:>6.3g}, Ei(x) = {Ei(x):>10.3g} : it goes to 0 quite fast!")
for x in np.linspace(1, 800, 9):
print(f"For x = {x:>6.3g}, Ei(x) = {Ei(x):>10.3g} : it goes to +oo quite fast!")
We can check that $\mathrm{Ei}(x)\to-\infty$ for $x\to0^-$ and $x\to0^+$:
for x in -1/np.logspace(1, 20, 10):
print(f"For x = {x:>10.3g} --> 0^-, Ei(x) = {Ei(x):>5.3g} : it doesn't go to -oo numerically!")
for x in 1/np.logspace(1, 20, 10):
print(f"For x = {x:>8.3g} --> 0^+, Ei(x) = {Ei(x):>5.3g} : it doesn't go to -oo numerically!")
Plots¶
And we can plot the $Ei(x)$ function, from $-1$ to $1$.
plt.plot(X, Y, 'b')
plt.title("The function $Ei(x)$ on $[-1,1]$")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()
Checking some inequalities¶
Let's check that $\forall x\in\mathbb{R}, \mathrm{Ei}(x) \leq \mathrm{e}^x$:
np.alltrue(Y <= np.exp(X))
We can check a tighter inequality, $\forall x\in\mathbb{R}, \mathrm{Ei}(x) \leq \mathrm{Ei}(-1) + (\mathrm{e}^x - \mathrm{e}) + (\mathrm{e} - \frac{1}{\mathrm{e}})$.
It is indeed tighter, as the constant on the right-hand side is non-negative:
Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1))
upper_bound = np.exp(X) + (Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1)))
np.alltrue(Y <= upper_bound)
plt.plot(X, Y, 'b')
plt.plot(X, np.exp(X), 'r--')
plt.plot(X, np.exp(X) + (Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1))), 'g--')
plt.title("The function $Ei(x)$ and upper-bound $e^x$ and $e^x + Ei(-1) - 1/e$")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()
We can check a tighter inequality, $\forall t\geq1, \forall x\geq1, \mathrm{Ei}(x) \leq \mathrm{Ei}(t) + \frac{\mathrm{e}^x - \mathrm{e}^{t}}{t}$.
e = np.exp(1)
upper_bound_cst = lambda t: Ei(t) - np.exp(t)/t
upper_bound_t = lambda t, X: Ei(t) + (np.exp(X) - np.exp(t))/t
upper_bound_cst(1)
upper_bound_cst(e)
upper_bound_cst(2*e)
X_4 = np.linspace(1, 2*e, 1000)
Y_4 = Ei(X_4)
def check_upper_bound(t):
upper_bound_4 = upper_bound_t(t, X_4)
return np.alltrue(Y_4 <= upper_bound_4)
check_upper_bound(1)
check_upper_bound(e)
check_upper_bound(2*e)
def see_upper_bound(t, xmax, onlylast=False):
X_4 = np.linspace(1, xmax, 1000)
Y_4 = Ei(X_4)
plt.plot(X_4, Y_4, 'b', label='Ei(x)')
upper_bound_4 = upper_bound_t(t, X_4)
plt.plot(X_4, upper_bound_4, 'y--', label='$Ei(t) + (e^x - e^t)/t$ for t = %.3g' % t)
if not onlylast:
plt.plot(X_4, np.exp(X_4), 'r--', label='$e^x$')
plt.plot(X_4, np.exp(X_4) + (Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1))), 'g--', label='$e^x + Ei(-1) - 1/e$')
plt.title("The function $Ei(x)$ and upper-bounds $e^x$ and $e^x + Ei(-1) - 1/e$ and $Ei(t) + (e^x - e^t)/t$ for t = %.3g" % t)
else:
plt.title("The function $Ei(x)$ and upper-bound $Ei(t) + (e^x - e^t)/t$ for t = %.3g" % t)
plt.legend()
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()
t = 1
see_upper_bound(t, 2*e)
t = 2
see_upper_bound(t, 2*e)
t = e
see_upper_bound(t, 2*e)
t = 2*e
see_upper_bound(t, t, onlylast=True)
t = 3*e
see_upper_bound(t, t, onlylast=True)
t = 4*e
see_upper_bound(t, t, onlylast=True)
I = lambda t: Ei(t) - Ei(-t)
I(1)
e - 1/e
assert I(1) < e - 1/e
I(e)
Other plots¶
X = np.linspace(1e-3, 2*e, 1000) # 1000 points
Y = Ei(X)
plt.plot(X, Y, 'b')
plt.title("The function $Ei(x)$ on $[0, e^2]$")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()
