How to calculate $\exp(-x)$ using Taylor series

Question

We know that the Taylor series expansion of $e^x$ is \begin{equation} e^x = \sum\limits_{i=1}^{\infty}\frac{x^{i-1}}{(i-1)!}. \end{equation}

If I have to use this formula to evaluate $e^{-20}$, how should I check for convergence?

MATLAB gives $e^{-20} = 2.0612e-009$. But when I use this formula to calculate it, I get $4.1736e-009$. My answer becomes more inaccurate as I take $x = -30, -40 \cdots$. The reason is obvious. This is because as the number of terms in the Taylor series increases, the factorial becomes more and more inaccurate using MATLAB.

1. How can I use Taylor series to get accurate values of $e^x$ in MATLAB?
2. How does MATLAB built in command "exp()" calculate the exponential?

Answer 1

The Problem

The main problem when computing $e^{-20}$ is that the terms of the series grow to $\frac{20^{20}}{20!}\approx43099804$ before getting smaller. Then the sum must cancel to be $\approx2.0611536\times10^{-9}$. In a floating point environment, this means that $16$ digits of accuracy in the sum are being thrown away due to the precision of the large numbers. This is the number of digits of accuracy of a double-precision floating point number ($53$ bits $\sim15.9$ digits).

For example, the RMS error in rounding $\frac{20^{20}}{20!}$, using double precision floating point arithmetic, would be $\sim\frac{20^{20}}{20!}\cdot2^{-53}/\sqrt3\approx3\times10^{-9}$. Since the final answer is $\approx2\times10^{-9}$, we lose all significance in the final answer simply by rounding that one term in the sum.

The problem gets worse with larger exponents. For $e^{-30}$, the terms grow to $\frac{30^{30}}{30!}\approx776207020880$ before getting smaller. Then the sum must cancel to be $\approx9.35762296884\times10^{-14}$. Here we lose $25$ digits of accuracy. For $e^{-40}$, we lose $33$ digits of accuracy.

---

A Solution

The usual solution is to compute $e^x$ and then use $e^{-x}=1/e^x$. When computing $e^x$, the final sum of the series is close in precision to the largest term of the series. Very little accuracy is lost.

For example, the RMS error in computing $e^{20}$ or $e^{-20}$, using double precision floating point arithmetic, would be $\sim8\times10^{-9}$; the errors are the same because both sums use the same terms, just with different signs. However, this means that using Taylor series, $$ \begin{align} e^{20\hphantom{-}}&=4.85165195409790278\times10^8\pm8\times10^{-9}\\ e^{-20}&=2\times10^{-9}\pm8\times10^{-9} \end{align} $$ Note that the computation of $e^{-20}$ is completely insignificant. On the other hand, taking the reciprocal of $e^{20}$, we get $$ e^{-20}=2.061153622438557828\times10^{-9}\pm3.4\times10^{-26} $$ which has almost $17$ digits of significance.

Answer 2

1) The regular double precision floating point arithmetic of Matlab is not sufficient to precisely calculate partial sums of this power series. To overcome that limitation, you can use the exact symbolic computation capabilities of the Symbolic Math Toolbox. This code

x = sym(-20);
i = sym(1 : 100);
expx = sum(x .^ (i - 1) ./ factorial(i - 1))

sums the first 100 terms of the series, and yields the exact (for the partial sum) result $$ \tiny\frac{20366871962240780739193874417376755657912596966525153814418798643652163252126492695723696663600640716177}{9881297415446727147594496649775206852319571477668037853762810667968023095834839075329261976769165978884198811117} $$ whose numeric counterpart (computed by double()) is 2.06115362243856e-09.

2) The Matlab function exp() is most likely a wrapper for C code, which presumably references the standard C library function exp().

You can have a look at the code of the widely used GNU C Library math functions, where you'll find that most functions are themselves wrappers around machine-specific implementations. For the i386 architecture, the path leads ultimately to the function __ieee754_exp which implements the exponential in inline assembler code (but don't ask me how to decipher this). On other architectures, there might be a single machine instruction that does the job.

"Most likely", "probably", "presumably" because Matlab is closed source and therefore definitive knowledge only exists within The Mathworks.

Answer 3

A standard technique, not particular to Matlab, is to exploit a function's symmetries to map its inputs to a more manageable range. The exponential function satisfies the symmetry $e^{a+b}=e^ae^b$, and it's easy to multiply floating-point numbers by powers of two. So do this: $$\begin{eqnarray} e^{-20}&=&2^{-20\lg e}\\ &\approx&2^{-28.854}\\ &=&2^{-29}\times2^{0.146}\\ &=&2^{-29}\times e^{0.146\log2}\\ &\approx&2^{-29}\times e^{0.101} \end{eqnarray}$$ Now you can easily evaluate $e^{0.101\ldots}$ using the Taylor series.