# coding: utf-8 # # Table of Contents #
# # Python implementation of [the Exponential Integral](https://en.wikipedia.org/wiki/Exponential_integral) function # # This small notebook is a [Python 3](https://www.python.org/) 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. $$ # In[25]: import numpy as np import matplotlib.pyplot as plt # In[26]: import seaborn as sns sns.set(context="notebook", style="darkgrid", palette="hls", font="sans-serif", font_scale=1.4) # In[27]: 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`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function. # In[28]: from scipy.integrate import quad # need only 1 function # First definition is the simplest: # In[29]: @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: # In[30]: @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: # In[31]: from numpy.linalg import norm # In[32]: X = np.linspace(-1, 1, 1000) # 1000 points Y = Ei(X) Y_2 = Ei_2(X) # In[33]: 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: # In[34]: get_ipython().run_line_magic('timeit', 'Y = Ei(X)') get_ipython().run_line_magic('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`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.expi.html#scipy.special.expi) # # The $\mathrm{Ei}$ function is also implemented as [`scipy.special.expi`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.expi.html#scipy.special.expi): # In[35]: from scipy.special import expi # In[36]: Y_3 = expi(X) # In[37]: np.allclose(Y, Y_3) # The difference is not too large: # In[38]: np.max(np.abs(Y - Y_3)) # In[39]: 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$. # In[40]: Ei(1) # In[41]: from scipy.optimize import root # In[42]: res = root(Ei, x0=1) res # In[43]: 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$: # In[44]: 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!") # In[45]: 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^+$: # In[46]: 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!") # In[47]: 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$. # In[48]: 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$: # In[49]: 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: # In[59]: Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1)) # In[60]: upper_bound = np.exp(X) + (Ei(-1) + (-np.exp(1)) + (np.exp(1) - np.exp(-1))) np.alltrue(Y <= upper_bound) # In[61]: 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}$. # In[76]: 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) # In[91]: 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) # In[107]: 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() # In[108]: t = 1 see_upper_bound(t, 2*e) # In[101]: t = 2 see_upper_bound(t, 2*e) # In[95]: t = e see_upper_bound(t, 2*e) # In[109]: t = 2*e see_upper_bound(t, t, onlylast=True) # In[110]: t = 3*e see_upper_bound(t, t, onlylast=True) # In[111]: t = 4*e see_upper_bound(t, t, onlylast=True) # In[113]: I = lambda t: Ei(t) - Ei(-t) I(1) e - 1/e assert I(1) < e - 1/e I(e) # ## Other plots # In[98]: 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() # ## Conclusion # # That's it, see [this page](https://en.wikipedia.org/wiki/Exponential_integral) or [this one](http://mathworld.wolfram.com/ExponentialIntegral.html) for more details on this function $\mathrm{Ei}(x)$.