{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Python implementation of [the Exponential Integral](https://en.wikipedia.org/wiki/Exponential_integral) function\n", "\n", "This small notebook is a [Python 3](https://www.python.org/) implementation of the Exponential Integral function, $Ei(x)$, defined like this:\n", "\n", "$$ \\forall x\\in\\mathbb{R}\\setminus\\{0\\},\\;\\; \\mathrm{Ei}(x) := \\int_{-\\infty}^x \\frac{\\mathrm{e}^u}{u} \\; \\mathrm{d}u. $$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import seaborn as sns\n", "sns.set(context=\"notebook\", style=\"darkgrid\", palette=\"hls\", font=\"sans-serif\", font_scale=1.4)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "import matplotlib as mpl\n", "mpl.rcParams['figure.figsize'] = (19.80, 10.80)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two implementations\n", "As one can show mathematically, there is another equivalent definition (the one used on Wikipedia):\n", "\n", "$$ \\forall x\\in\\mathbb{R}\\setminus\\{0\\},\\;\\; \\mathrm{Ei}(x) := - \\int_{-x}^{\\infty} \\frac{\\mathrm{e}^{-t}}{t} \\; \\mathrm{d}t. $$\n", "\n", "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$.\n", "\n", "We use the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import quad # need only 1 function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First definition is the simplest:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "@np.vectorize\n", "def Ei(x, minfloat=1e-7, maxfloat=10000):\n", " \"\"\"Ei integral function.\"\"\"\n", " minfloat = min(np.abs(x), minfloat)\n", " maxfloat = max(np.abs(x), maxfloat)\n", " def f(t):\n", " return np.exp(t) / t\n", " if x > 0:\n", " return (quad(f, -maxfloat, -minfloat)[0] + quad(f, minfloat, x)[0])\n", " else:\n", " return quad(f, -maxfloat, x)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other definition is very similar:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "@np.vectorize\n", "def Ei_2(x, minfloat=1e-7, maxfloat=10000):\n", " \"\"\"Ei integral function.\"\"\"\n", " minfloat = min(np.abs(x), minfloat)\n", " maxfloat = max(np.abs(x), maxfloat)\n", " def f(t):\n", " return np.exp(-t) / t\n", " if x > 0:\n", " return -1.0 * (quad(f, -x, -minfloat)[0] + quad(f, minfloat, maxfloat)[0])\n", " else:\n", " return -1.0 * quad(f, -x, maxfloat)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking the two versions\n", "We can quickly check that the two are equal:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "from numpy.linalg import norm" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "X = np.linspace(-1, 1, 1000) # 1000 points\n", "Y = Ei(X)\n", "Y_2 = Ei_2(X)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Two versions of Ei(x) are indeed equal for 1000 values.\n" ] } ], "source": [ "assert np.allclose(Y, Y_2)\n", "print(f\"Two versions of Ei(x) are indeed equal for {len(X)} values.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compare which is fastest to evaluate:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.72 s ± 176 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "1.62 s ± 18.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit Y = Ei(X)\n", "%timeit Y_2 = Ei_2(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "They both take about the same time, but the second implementation seems (slightly) faster." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparison with [`scipy.special.expi`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.expi.html#scipy.special.expi)\n", "\n", "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):" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "from scipy.special import expi" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "Y_3 = expi(X)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(Y, Y_3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The difference is not too large:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.0000000611197777e-07" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.max(np.abs(Y - Y_3))" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our version of Ei(x) is the same as the one in scipy.special.expi (1000 values).\n" ] } ], "source": [ "assert np.allclose(Y, Y_3, rtol=1e-6, atol=1e-6)\n", "print(f\"Our version of Ei(x) is the same as the one in scipy.special.expi ({len(X)} values).\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Special values\n", "We can compute some special values, like $\\mathrm{Ei}(1)$ and solving (numerically) $\\mathrm{Ei}(x)=0$." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(1.89511762)" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Ei(1)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import root" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " fjac: array([[-1.]])\n", " fun: array([-3.55271368e-15])\n", " message: 'The solution converged.'\n", " nfev: 9\n", " qtf: array([-2.78710388e-12])\n", " r: array([-3.89621516])\n", " status: 1\n", " success: True\n", " x: array([0.37250746])" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res = root(Ei, x0=1)\n", "res" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The approximate solution to Ei(x)=0 is x0 = 0.37250746211322744 (for which Ei(x)=[-3.55271368e-15])...\n" ] } ], "source": [ "print(f\"The approximate solution to Ei(x)=0 is x0 = {res.x[0]} (for which Ei(x)={res.fun})...\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Limits" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check that $\\mathrm{Ei}(x)\\to0$ for $x\\to-\\infty$ and $\\mathrm{Ei}(x)\\to+\\infty$ for $x\\to\\infty$:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For x = -1, Ei(x) = -0.219 : it goes to 0 quite fast!\n", "For x = -112, Ei(x) = -1.17e-54 : it goes to 0 quite fast!\n", "For x = -223, Ei(x) = -4.27e-103 : it goes to 0 quite fast!\n", "For x = -334, Ei(x) = -2e-151 : it goes to 0 quite fast!\n", "For x = -445, Ei(x) = -1.05e-199 : it goes to 0 quite fast!\n", "For x = -556, Ei(x) = -5.85e-248 : it goes to 0 quite fast!\n", "For x = -667, Ei(x) = -3.39e-296 : it goes to 0 quite fast!\n", "For x = -778, Ei(x) = -0 : it goes to 0 quite fast!\n", "For x = -889, Ei(x) = -0 : it goes to 0 quite fast!\n", "For x = -1e+03, Ei(x) = -0 : it goes to 0 quite fast!\n" ] } ], "source": [ "for x in -np.linspace(1, 1000, 10):\n", " print(f\"For x = {x:>6.3g}, Ei(x) = {Ei(x):>10.3g} : it goes to 0 quite fast!\")" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For x = 1, Ei(x) = 1.9 : it goes to +oo quite fast!\n", "For x = 101, Ei(x) = 6.46e+41 : it goes to +oo quite fast!\n", "For x = 201, Ei(x) = 7.66e+84 : it goes to +oo quite fast!\n", "For x = 301, Ei(x) = 1.21e+128 : it goes to +oo quite fast!\n", "For x = 400, Ei(x) = 2.15e+171 : it goes to +oo quite fast!\n", "For x = 500, Ei(x) = 4.09e+214 : it goes to +oo quite fast!\n", "For x = 600, Ei(x) = 8.08e+257 : it goes to +oo quite fast!\n", "For x = 700, Ei(x) = 1.64e+301 : it goes to +oo quite fast!\n", "For x = 800, Ei(x) = inf : it goes to +oo quite fast!\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:7: RuntimeWarning: overflow encountered in exp\n", " import sys\n", "/usr/local/lib/python3.6/dist-packages/numpy/lib/function_base.py:2831: RuntimeWarning: invalid value encountered in ? (vectorized)\n", " outputs = ufunc(*inputs)\n" ] } ], "source": [ "for x in np.linspace(1, 800, 9):\n", " print(f\"For x = {x:>6.3g}, Ei(x) = {Ei(x):>10.3g} : it goes to +oo quite fast!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check that $\\mathrm{Ei}(x)\\to-\\infty$ for $x\\to0^-$ and $x\\to0^+$:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For x = -0.1 --> 0^-, Ei(x) = -1.82 : it doesn't go to -oo numerically!\n", "For x = -0.000774 --> 0^-, Ei(x) = -6.59 : it doesn't go to -oo numerically!\n", "For x = -5.99e-06 --> 0^-, Ei(x) = -11.4 : it doesn't go to -oo numerically!\n", "For x = -4.64e-08 --> 0^-, Ei(x) = -16.3 : it doesn't go to -oo numerically!\n", "For x = -3.59e-10 --> 0^-, Ei(x) = -21.2 : it doesn't go to -oo numerically!\n", "For x = -2.78e-12 --> 0^-, Ei(x) = -26 : it doesn't go to -oo numerically!\n", "For x = -2.15e-14 --> 0^-, Ei(x) = -30.8 : it doesn't go to -oo numerically!\n", "For x = -1.67e-16 --> 0^-, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n", "For x = -1.29e-18 --> 0^-, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n", "For x = -1e-20 --> 0^-, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n", " If increasing the limit yields no improvement it is advised to analyze \n", " the integrand in order to determine the difficulties. If the position of a \n", " local difficulty can be determined (singularity, discontinuity) one will \n", " probably gain from splitting up the interval and calling the integrator \n", " on the subranges. Perhaps a special-purpose integrator should be used.\n", " warnings.warn(msg, IntegrationWarning)\n" ] } ], "source": [ "for x in -1/np.logspace(1, 20, 10):\n", " print(f\"For x = {x:>10.3g} --> 0^-, Ei(x) = {Ei(x):>5.3g} : it doesn't go to -oo numerically!\")" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For x = 0.1 --> 0^+, Ei(x) = -1.62 : it doesn't go to -oo numerically!\n", "For x = 0.000774 --> 0^+, Ei(x) = -6.59 : it doesn't go to -oo numerically!\n", "For x = 5.99e-06 --> 0^+, Ei(x) = -11.4 : it doesn't go to -oo numerically!\n", "For x = 4.64e-08 --> 0^+, Ei(x) = -16.3 : it doesn't go to -oo numerically!\n", "For x = 3.59e-10 --> 0^+, Ei(x) = -21.2 : it doesn't go to -oo numerically!\n", "For x = 2.78e-12 --> 0^+, Ei(x) = -26 : it doesn't go to -oo numerically!\n", "For x = 2.15e-14 --> 0^+, Ei(x) = -30.8 : it doesn't go to -oo numerically!\n", "For x = 1.67e-16 --> 0^+, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n", "For x = 1.29e-18 --> 0^+, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n", "For x = 1e-20 --> 0^+, Ei(x) = -31.9 : it doesn't go to -oo numerically!\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n", " If increasing the limit yields no improvement it is advised to analyze \n", " the integrand in order to determine the difficulties. If the position of a \n", " local difficulty can be determined (singularity, discontinuity) one will \n", " probably gain from splitting up the interval and calling the integrator \n", " on the subranges. Perhaps a special-purpose integrator should be used.\n", " warnings.warn(msg, IntegrationWarning)\n" ] } ], "source": [ "for x in 1/np.logspace(1, 20, 10):\n", " print(f\"For x = {x:>8.3g} --> 0^+, Ei(x) = {Ei(x):>5.3g} : it doesn't go to -oo numerically!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plots\n", "And we can plot the $Ei(x)$ function, from $-1$ to $1$." ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[