Numerical sensitivity + floating point roundoff errors for a 2013 computational physics course

Gauss.md

Suppose little Gauss lived in the modern age. Little Gauss’ teacher wanted to surf the internet, so he assigned all of his students the following integral to evaluate:

$$ \int_0^1 x^{100} e^{x-1} dx $$

Being the clever alter-ego of the boy who immediately saw $\sum^n_k k = {n+1 \choose 2}$, little Gauss constructed a sequence

$$ s_k = \int_0^1 x^k e^{x-1} dx $$

and with a little bit of clever manipulation (integration by parts), he found that, using $u_k = x^k$, $dv = e^{x-1}dx$, $v = e^{x-1}$

$$ \begin{aligned} s_k &= \int_0^1 u_kdv = \left.u_kv\right|_0^1 - \int_0^1 vdu_k \\ &= \left. x^ke^{x-1}\right|_0^1 - k\int_0^1x^{k-1}e^{x-1}dx \\ &= \left. x^ke^{x-1}\right|_0^1 - ks_{k-1} \\ &= 1 - ks_{k-1} \end{aligned} $$

and

$$ s_0 = \int_0^1 e^{x-1}dx = 1 - 1/e $$

Why, this is a simple recursive function! Little Gauss was absolutely thrilled, he has at his disposal a programmable calculator capable of python (because he’s Gauss, he can have whatever the fuck he wants), and he quickly coded up the recurrence:

import math
def s(k):
    if k == 0:
        return 1 - math.exp(-1)
    else:
        return 1 - k*s(k-1)

He verified the first few cases by hand, and upon finding his code satisfactory, he computed the 100th element of the sequence as -1.1599285429663592e+141 and turned in the first few digits.

His teacher glanced at his solution, and knowing that there’s no way little Gauss could have done his school work with such proficiency, immediately declared it wrong. Upon repeated appeal, the teach finally relented and looked up the solution in his solution manual and, bewildered… again told little Gauss that he was WAAAAY off. Unfortunately for our tragic hero, this would not be a modern day retelling of a clever little boy who outsmarted his teacher. No, this is a tragic story of a clever little boy who succumbed to a fatal case of the roundoff bugs.

---

Suppose $f(x) = x^{100}e^{x-1}$, because $0 \le e^{x-1} \le 1$ inside $0 \le x \le 1$, then $\frac{x^{100}}{e} \le x^{100}e^{x-1} \le x^{100}$, so

$$ \int_0^1 \frac{x^{100}}{e} dx \le \int_0^1 x^{100}e^{x-1} dx \le \int_0^1 x^{100} dx $$

or

$$ \frac{1}{e(100 + 1)} \le s_{100} \le \frac{1}{100 + 1} $$

and in general

$$ \frac{1}{e(k+1)} \le s_k \le \frac{1}{k+1} $$

of course $s_{100}$ can’t be on the order of $-10^{141}$

So what went wrong? Well, whenever you’re in trouble, just make a plot.

Source Source

Everything seems to be going fine until around $k=17$.

It turns out that there was nothing wrong with little Gauss’ method and the integral is perfectly well-behaved. The culprit lies in the fact that $s_0$ can only be represented approximately in the computer. As we know already, $s_0 = 1 - e^{-1}$ can only be represented approximately. We’re actually giving using an initial $\hat s_0 = s_0 + \epsilon$ for some tiny perturbation $\epsilon$. Now, let’s see what the computer is computing once we unroll the recursion (we’ll use the notation $\hat s_k$ to be the $s_k$ term plus the propagated floating point error):

$$ \begin{aligned} \hat s_0 &= s_0 + \epsilon \\ \hat s_1 &= 1 - \hat s_0 = s_1 - \epsilon \\ \hat s_2 &= 1 - 2\hat s_1 = s_2 + 2! \epsilon \\ \hat s_3 &= 1 - 3\hat s_2 = s_3 - 3! \epsilon \\ &\vdots\\ \hat s_k &= s_{k} \pm k!\epsilon \end{aligned} $$

The tiny perturbation in the computed value of $\hat s_0$ propagated exponentially throughout the computation and at the $k^{th}$ step, manifests itself as $k! \epsilon$. Even at $k=19$, we will see around $10^{17}\epsilon$ added into the calculation. For $k=100$, the answer will have an additional factor of $10^{158}\epsilon$. No matter how small your initial $\epsilon$ is, $100! \epsilon$ will definitely be noticeable.

Little Gauss definitely should have paid attention to the lessons on round-off errors.

Problem: How would you turn this into a stable algorithm?