Exponential integral representation

Question

According to exponential integral eqn. (8) $\; E_{1}(x) \;$ can be represented by:

$$ E_1(x)= - \gamma - \ln(x) - \sum _{n=1}^{\infty } \frac{(-1)^n x^n}{n n!} $$

where $\gamma$ is the Euler-Mascheroni constant.

I observed that this representation is valid or give accurate results for lower values of $x$ only, i.e. $x=0.001, ...,20$. However, for large values of $x$ this representation gives wrong result.

For example, for $x=0.1$, the result is $1.82292$, which is correct. For $x=25$, this representation gives $-0.00002210807401$, which is not correct, it should be $5.348899755*10^{-13}$.

Any explanation of this?

Answer 1

Generally, power series representation is not a good choice for larger input, even if the radius of convergence if infinite.

In this case, you need about 90 leading terms just simply to achieve the desired order of magnitude (and more terms for the accurate digits). Indeed, Mathematica 10.2 confirms that if we let

$$ S_n = -\gamma - \log 25 - \sum_{k=1}^{n} \frac{(-1)^k 25^k}{k!k}, $$

then you can see that we need about 120 leading terms to achieve 20-digit accuracy:

\begin{align*} S_{10} &\approx -1.8036638952044776218 \times 10^{6} \\ S_{20} &\approx -1.0045038653588617198 \times 10^{8} \\ S_{30} &\approx -4.8147011301872907116 \times 10^{7} \\ S_{40} &\approx -9.5062897788081758876 \times 10^{5} \\ S_{50} &\approx -1.6908961484269455556 \times 10^{3} \\ S_{60} &\approx -4.3435091996295920753 \times 10^{-1} \\ S_{70} &\approx -2.2108074007813609212 \times 10^{-5} \\ S_{80} &\approx -2.7925545196351250856 \times 10^{-10} \\ S_{90} &\approx 5.3384543545972366461 \times 10^{-13} \\ S_{100} &\approx 5.3488997421949443337 \times 10^{-13} \\ S_{110} &\approx 5.3488997553402104337 \times 10^{-13} \\ S_{110} &\approx 5.3488997553402166403 \times 10^{-13} \end{align*}

This also says that we have horrible cancellations between leading terms, which are much larger compared to the magnitude of $E_1(25)$. So if you use usual floating type variables, then you will certainly lose all the significant digits. For example, double-precision floating type gives 15-17 significant digits, which is clearly inappropriate in our case. I suspect this is the reason why you get different values.

For large $x$, as @gammatester pointed out, other methods are more adequate for numerical purpose.

Answer 2

Without seeing your code it is difficult. But I guess it is catastrophic cancellation. $E_1(25) \approx 5.3488997553\times 10^{-13}$ and the single terms are much larger, e.g. the for $n=10$ the value is $2628070.75729$ and for $n=24$ the term is $238585088.1445781!$ You will get similar problems is you compute $e^{-x}$ or $\sin(x)$ for $x=25$ using the Taylor series.

If you are interesting in actually computing $E_1$ (and not only to get insight in the described problem) you can use the continued fraction http://dlmf.nist.gov/6.9 for say $x>1$.

Answer 3

Looking here, there is an interesting expansion for large arguments $$E_1(x)=\frac{e^{-x}}{x}\sum_{k=0}^\infty \frac {k!}{(-x)^k}$$ For $x=25$, using $n$ terms, we have $$S_1=5.332970444146184 \times 10^{-13}$$ $$S_2=5.350747012293338 \times 10^{-13}$$ $$S_3=5.348613824115680 \times 10^{-13}$$ $$S_4=5.348955134224105 \times 10^{-13}$$ $$S_5=5.348886872202420\times 10^{-13}$$ $$S_6=5.348903255087624\times 10^{-13}$$ $$S_7=5.348898667879767\times 10^{-13}$$ $$S_8=5.348900135786281\times 10^{-13}$$ $$S_9=5.348899607339937\times 10^{-13}$$

Edit

As @tired commented, we must take care about the fact that this asymptotic expansion diverges for large values of $k$. In the present case, we start with trouble around $k=50$.