<a href="https://colab.research.google.com/github/stephenbeckr/numerical-analysis-class/blob/master/Demos/Ch1_ExponentialSummation.ipynb" target="_parent"></a>

In [Ch1_SymbolicTaylorSeries](Ch1_SymbolicTaylorSeries.ipynb) we saw how to use the Taylor remainder theorem to estimate error in an expansion, and this is useful if you want to calculate a function like sine or cosine using only additions/subtractions and multiplications/division.

Let's try that with the exponential function.  This is a classic example in numerical analysis


```
# Load the package
import sympy as sym
from sympy import init_printing
init_printing()  # This will make output look very nice!

import math  # has things like math.pi (which is the numerical value of pi)
```


```
# First, all symbolic variables must be declared. This is in contrast to Mathematica
x = sym.Symbol('x')

# Recall the Taylor series for exp(x)... very simple!
sym.series( sym.exp(-x), x )
```

Let's try to find $e^{-30}$ with about 14 digits of precision.  Recally that Taylor's remainder theorem says that the error in the Taylor series (with $n-1$ terms) is the $n$th term but evaluate the derivative at some $\xi \in [x_0,x]$.  We'll expand around $0$, so $x_0=0$ (we do this because we *know* that $e^0=1$; we could choose a different $x_0$ to expand about if we knew what $e^{x_0}$ is).

For this remainder term, the $n$th derivative of $e^{-x}$ is $(-1)^ne^{-x}$, and since our target $x=30$ is positive, we know $\xi \ge 0$ and hence $e^{-\xi} \le 1$, so we can bound this derivative term by $1$ (in absolute value).  We also need to include the coefficient $1/n!$ and the polynomial term $(x-x_0)^n$.


```
X   = -30

# Taylor's remainder theorem says that the remainder term at the (n-1) step is the nth Taylor series term but evaluated at some \xi
Rem = lambda n : abs( X**n / math.factorial(n) )
for n in range(15):
  print( Rem(n) )
  # When n is large enough, eventually this sequence will start to decrease
```

    1.0
    30.0
    450.0
    4500.0
    33750.0
    202500.0
    1012500.0
    4339285.714285715
    16272321.42857143
    54241071.428571425
    162723214.2857143
    443790584.41558444
    1109476461.038961
    2560330294.7052946
    5486422060.082774


So if we want 14 digits, let's calculate a bound on the number of Taylor series we need:


```
tol = 1.0e-14
n = 0
while Rem(n) > tol and n < 1e8 : # this last check is to prevent infinite loops, which are a possibility due to bugs
  n += 1
N = n-1
print( "We need at least", N, "terms, since Rem(n+1) is %.2e" % Rem(N+1) )
```

    We need at least 106 terms, since Rem(n+1) is 9.19e-15


Now, let's sum the Taylor series. We should have an error of 1e-14, right?


```
e = 0.
for n in range(N+1):
  e += X**n / math.factorial(n)
  if ( not (n % 10 ) ) or ( n > N-10 ):  # save some output...
    print("Iteration %3d"% n, " estimate is %29.15f"% e, "and error is %9.3e" % abs(e - math.exp(X)) )
```

    Iteration   0  estimate is             1.000000000000000 and error is 1.000e+00
    Iteration  10  estimate is     121254849.571428596973419 and error is 1.213e+08
    Iteration  20  estimate is   85291712210.032180786132812 and error is 8.529e+10
    Iteration  30  estimate is  384842612535.907775878906250 and error is 3.848e+11
    Iteration  40  estimate is   63336540117.266479492187500 and error is 6.334e+10
    Iteration  50  estimate is     878229229.277898311614990 and error is 8.782e+08
    Iteration  60  estimate is       1685584.300283375196159 and error is 1.686e+06
    Iteration  70  estimate is           622.524572911353516 and error is 6.225e+02
    Iteration  80  estimate is             0.055867710574729 and error is 5.587e-02
    Iteration  90  estimate is            -0.000084070677839 and error is 8.407e-05
    Iteration  97  estimate is            -0.000085530630269 and error is 8.553e-05
    Iteration  98  estimate is            -0.000085530022814 and error is 8.553e-05
    Iteration  99  estimate is            -0.000085530206891 and error is 8.553e-05
    Iteration 100  estimate is            -0.000085530151668 and error is 8.553e-05
    Iteration 101  estimate is            -0.000085530168071 and error is 8.553e-05
    Iteration 102  estimate is            -0.000085530163247 and error is 8.553e-05
    Iteration 103  estimate is            -0.000085530164652 and error is 8.553e-05
    Iteration 104  estimate is            -0.000085530164246 and error is 8.553e-05
    Iteration 105  estimate is            -0.000085530164362 and error is 8.553e-05
    Iteration 106  estimate is            -0.000085530164330 and error is 8.553e-05


... but our error is in fact much worse. What went wrong?  Our estimate of the error was due to our mathematical approximation of the exponential, but we didn't include numerical error due to the summation.  We are adding and subtracting a lot of numbers of similar sizes, so we loose precision.

Attempted fix #1: use better summation strategies.  Unfortunately, we'll see that this doesn't quite fix it for us.
- Fix 1a: use Python's builtin `fsum`
- Fix 1b: use [Kahan summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm), which is a trick that effectively let's us work with twice as many digits, or use another good summing algorithm. 

---




```
# Fix 1a.
# See https://www.math.ubc.ca/~pwalls/math-python/python/sequences/
series = [ X**n / math.factorial(n) for n in range(N+1) ]  # "lazy", or "iterable". This is a list comprehension
e = math.fsum( series )
abs( e - math.exp(X) )
```


```
# Kahan Summation, though doesn't work (Python might be optimizing?)
e = 0.
accum = 0.
for n in range(N+1):
  a = X**n / math.factorial(n) - accum
  b = e + a
  accum = (b-e) -a
  e = b
  if n > N-4:
    print("Iteration %2d"% n, " estimate is %21.15f"% e, "and error is %9.3e" % abs(e - math.exp(X)) )
```

    Iteration 103  estimate is    -0.000085530164652 and error is 8.553e-05
    Iteration 104  estimate is    -0.000085530164246 and error is 8.553e-05
    Iteration 105  estimate is    -0.000085530164362 and error is 8.553e-05
    Iteration 106  estimate is    -0.000085530164330 and error is 8.553e-05


Attempted fix #2: use some math.
- Fix 2a: split into positive and negative terms, add those separately, then combine at the end. Good idea, but again, doesn't quite fix it.
- Fix 2b: rewrite $e^{-x} = \frac{1}{e^x} = \frac{1}{1+x+x^2/2+\ldots}$.  Finally, we have a good solution!


```
# Fix 2a
seriesEven = [ X**n / math.factorial(n) for n in range(0,N+1,2) ] 
seriesOdd = [ X**n / math.factorial(n) for n in range(1,N+1,2) ] 
e = math.fsum( seriesEven ) + math.fsum( seriesOdd )
abs( e - math.exp(X) )
```


```
# Fix 2b
series = [ abs(X)**n / math.factorial(n) for n in range(N+1) ] 
e = 1./math.fsum( series )
abs( e - math.exp(X) )
```


```

```

fd5285bf0fd13b101cb925b9d1f7964169879c69
ipynb
Jupyter Notebook
Demos/Ch1_ExponentialSummation.ipynb
skhadem/numerical-analysis-class
a022fc2e73254800a1e193f94280446223ef01ae
BSD-3-Clause
2020-08-25T19:11:46.000Z
2021-11-12T21:49:44.000Z
Demos/Ch1_ExponentialSummation.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-09-01T21:44:12.000Z
2020-09-01T21:44:12.000Z
Demos/Ch1_ExponentialSummation.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-08-25T21:25:17.000Z
2021-03-03T04:14:26.000Z
Qwen/Qwen-72B
1. YES
2. YES