6
$\begingroup$

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?
$\endgroup$
  • 1
    $\begingroup$ If the ultimate goal is to approximate the exponential I would use a Pade' approximants (a rational approximation), combined with $e^{2x}=(e^x)^2$ to get values for larger values of $x$. $\endgroup$ – Pp.. Jan 28 '15 at 23:39
  • $\begingroup$ If you just want to see how the Taylor series behaves, then the error of an alternating series can be bounded by the absolute value of the next term of the series. $\endgroup$ – Pp.. Jan 28 '15 at 23:40
  • 1
    $\begingroup$ TS is good for small values of the argument. Otherwise you might end up having many extra terms. $\endgroup$ – Kaster Jan 28 '15 at 23:42
  • $\begingroup$ You may try summing backwards. $\endgroup$ – lhf Jan 29 '15 at 11:44
  • $\begingroup$ Oog. Why not start the summation at $0$? $\endgroup$ – Qiaochu Yuan Jan 30 '15 at 8:54
11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Would the downvoter care to comment? $\endgroup$ – robjohn Jan 30 '15 at 1:43
  • $\begingroup$ this is so nice an explanation!! When doing numerical calculation the formulas must be modified in a manner to avoid loss of significant digits and such that intermediate calculation does not require a considerably higher precision than that desired in the final answer. $\endgroup$ – Paramanand Singh Feb 19 '15 at 10:48
3
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You mention using "variable precision arithmetic" in Matlab but your example code is an exact symbolic calculation converted to floating-point. There appears to be no need for vpa for the ranges used. The OP also asked "how should I check for convergence?" which I think was the impetus for the question. $\endgroup$ – horchler Jan 29 '15 at 3:15
  • $\begingroup$ @horchler, correct with the vpa thing, I got the two confused, edited. – As for convergence, I think this was a minor part of the question made irrelevant by the fact that his double precision computation couldn't converge anyway. 4.1e-9 isn't not-converged, its just off. Feel free to add a better answer though, or even to edit mine. :-) $\endgroup$ – A. Donda Jan 29 '15 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.