1.2 Fehlerfortpflanzung

Frage: Die Eingaben $x,y$ sind schon fehlerbehaftet (Messung, Rundung, Vorrechnung). Wie groß ist dann der Fehler im Ergebnis $z=x+y$, $z=x\cdot y$ oder allgemein $z=f(x_1,\dots,x_n)$?

Absoluter und relativer Fehler (Tafel)

Sei $\tilde x$ eine Näherung für den exakten Wert $x$.

$$\Delta x=\tilde x-x\quad\text{(absoluter Fehler)},\qquad \varepsilon_x=\frac{\Delta x}{x}=\frac{\tilde x-x}{x}\quad\text{(relativer Fehler, }x\neq0\text{)}.$$

Damit ist $\tilde x=x+\Delta x=x(1+\varepsilon_x)$.

Der relative Fehler ist die natürliche Größe für Gleitkommazahlen: $|\varepsilon_x|\approx10^{-k}$ heißt „etwa $k$ gültige Dezimalstellen“. Die Rundung macht laut Rechenmodell genau $|\varepsilon|\le\text{eps}$.

Addition $z=x+y$

$\Delta z=\tilde z-z=(x+\Delta x)+(y+\Delta y)-(x+y)=\Delta x+\Delta y.$

Absolut addieren sich die Fehler einfach. Relativ:

$$\varepsilon_z=\frac{\Delta x+\Delta y}{x+y}=\frac{x}{x+y}\,\varepsilon_x+\frac{y}{x+y}\,\varepsilon_y.$$

Die Faktoren $\frac{x}{x+y}$ und $\frac{y}{x+y}$ heißen Verstärkungsfaktoren. Entscheidend ist, ob sie klein bleiben.

1. Fall: gleiche Vorzeichen

Für $x,y>0$ (oder beide $<0$) gilt $0<\frac{x}{x+y},\frac{y}{x+y}<1$ und beide summieren sich zu 1. $\varepsilon_z$ ist also ein gewichtetes Mittel von $\varepsilon_x$ und $\varepsilon_y$:

$|\varepsilon_z|\le\max\big(|\varepsilon_x|,|\varepsilon_y|\big).$

⇒ Die Addition von Zahlen gleichen Vorzeichens ist stabil: Der relative Fehler wird nicht größer.

2. Fall: verschiedene Vorzeichen (Auslöschung)

Für $x>0>y$ ist $|x+y|=\big||x|-|y|\big|$, also

$$|\varepsilon_z|\le\frac{|x|}{\big||x|-|y|\big|}\,|\varepsilon_x|+\frac{|y|}{\big||x|-|y|\big|}\,|\varepsilon_y|.$$

Gilt $|x|\approx|y|$, dann ist der Nenner winzig und die Verstärkungsfaktoren werden beliebig groß. ⇒ Die Subtraktion fast gleich großer Zahlen ist instabil. Das heißt Auslöschung: Die führenden, übereinstimmenden Ziffern heben sich auf, übrig bleiben die hinteren, fehlerbehafteten.

Beispiel (Folie §1/4) mit $N=8$ Stellen: 📄 Folie §1/4

$$\begin{aligned} x&=0.98765444\ldots,\quad \tilde x=0.98765444\\ y&=0.98765333\ldots,\quad \tilde y=0.98765333\\ \tilde x-\tilde y&=0.00000111=1.11\cdot10^{-6}. \end{aligned}$$

Eingaben: relativer Fehler $\approx10^{-8}$ (8 gültige Stellen). Ergebnis: nur noch 3 gültige Stellen ($1.11$), denn die Stellen dahinter sind unbekannt. Der relative Fehler ist um den Faktor $\frac{0.98}{1.11\cdot10^{-6}}\approx10^{6}$ gewachsen.

Merke

Auslöschung erzeugt keinen neuen Fehler. Die Subtraktion selbst ist sogar exakt. Sie macht den schon vorhandenen Fehler von $x$ und $y$ relativ sichtbar, weil das Ergebnis viel kleiner ist als die Eingaben.

Allgemeine Funktion: totales Differential (Tafel)

Sei $f:D\subseteq\mathbb R^n\to\mathbb R$ differenzierbar, $z=f(x_1,\dots,x_n)$. Taylor um $x=(x_1,\dots,x_n)$ bis zur ersten Ordnung:

$$\tilde z=f(x+\Delta x)=f(x)+\sum_{i=1}^n\frac{\partial f}{\partial x_i}(x)\,\Delta x_i+\mathcal O\big(\|\Delta x\|^2\big).$$

Terme höherer Ordnung werden vernachlässigt (die Fehler sind klein, ihre Produkte noch viel kleiner):

$$\boxed{\;\Delta z\approx\sum_{i=1}^n\frac{\partial f}{\partial x_i}\,\Delta x_i,\qquad \varepsilon_z\approx\sum_{i=1}^n\underbrace{\frac{x_i}{f(x)}\frac{\partial f}{\partial x_i}}_{\text{Verstärkungsfaktor}}\,\varepsilon_{x_i}\;}$$

Die Verstärkungsfaktoren sind genau die relativen Konditionszahlen aus Kapitel 0. Für $f(x,y)=x+y$ kommt die Formel von oben heraus.

Multiplikation $z=x\cdot y$

$\frac{\partial f}{\partial x}=y$, $\frac{\partial f}{\partial y}=x$:

$\Delta z\approx y\,\Delta x+x\,\Delta y,\qquad\varepsilon_z\approx\frac{xy}{xy}\varepsilon_x+\frac{xy}{xy}\varepsilon_y=\varepsilon_x+\varepsilon_y.$

⇒ Die relativen Fehler addieren sich, die Faktoren sind immer 1. Multiplikation (und analog Division: $\varepsilon_{x/y}\approx\varepsilon_x-\varepsilon_y$) ist stabil.

Logarithmus $z=\ln x$

$f'(x)=1/x$:

$\Delta z\approx\frac{\Delta x}{x}=\varepsilon_x,\qquad\varepsilon_z\approx\frac{\varepsilon_x}{\ln x}.$

⇒ Kritisch für $x\approx1$: Dort ist $\ln x\approx0$ und der Faktor $1/\ln x$ explodiert. Beispiel: $x=1.001$, $\ln x\approx10^{-3}$, also wird $\varepsilon_x$ um den Faktor $1000$ verstärkt (3 Stellen verloren).

Tabelle elementarer Funktionen (Folie §1/5)

📄 Folie §1/5

$f(x)$	$\varepsilon_z\approx$	kritisch, wenn …
$x^\alpha$	$\alpha\,\varepsilon_x$	$\lvert\alpha\rvert$ groß
$\sqrt x$	$\tfrac12\,\varepsilon_x$	nie (halbiert den Fehler)
$e^x$	$x\,\varepsilon_x$	$\lvert x\rvert$ groß
$\ln x$	$\varepsilon_x/\ln x$	$x\approx1$
$\sin x$	$\varepsilon_x\cdot x\cot x$	$x\approx k\pi$, $k\neq0$
$\cos x$	$-\varepsilon_x\cdot x\tan x$	$x\approx\frac\pi2+k\pi$
$\arctan x$	$\varepsilon_x\cdot\frac{x}{(1+x^2)\arctan x}$	nie (Faktor $\le1$)
$\tanh x$	$\varepsilon_x\cdot\frac{2x}{\sinh 2x}$	nie (Faktor $\le1$)

Ablesen geht immer gleich: Verstärkungsfaktor $=\dfrac{x\,f'(x)}{f(x)}$. Die Zeilen $x^\alpha$ und $e^x$ sind Ergänzung, der Rest steht so auf der Folie (dort für $\tanh$ als $\varepsilon_x\cdot x\cosh^{-2}x\coth x$, dasselbe).

Ergänzung: die Zeilen nachrechnen

Für $\sin$: $\frac{x\cos x}{\sin x}=x\cot x$. Bei $x\to0$ geht das gegen $1$ (harmlos), bei $x\to\pi$ gegen $\infty$, weil $\sin\pi=0$. Für $\tanh$: $\frac{x(1-\tanh^2x)}{\tanh x}=\frac{x}{\sinh x\cosh x}=\frac{2x}{\sinh2x}\le1$.

Zusätzliche Rundung im Ergebnis

Bisher wurde $z$ exakt aus den fehlerhaften Eingaben berechnet. Am Rechner wird das Ergebnis aber selbst noch gerundet (Rechenmodell): $\tilde z=\operatorname{rd}(\tilde x+\tilde y)=(\tilde x+\tilde y)(1+\tau_z)$, $|\tau_z|\le\text{eps}$. In erster Ordnung:

$$\begin{aligned} z=x+y:&\quad\Delta z\approx\Delta x+\Delta y+z\,\tau_z,\\ z=x\cdot y:&\quad\varepsilon_z\approx\varepsilon_x+\varepsilon_y+\tau_z,\\ z=f(x):&\quad\varepsilon_z\approx\frac{f'(x)\,\Delta x}{f(x)}+\tau_z. \end{aligned}$$

Jede Operation kostet also zusätzlich höchstens ein $\text{eps}$ relativ. Über einen ganzen Algorithmus mit $k$ Operationen wächst das im harmlosen Fall wie $k\cdot\text{eps}$. Gefährlich wird es nur, wenn irgendwo große Verstärkungsfaktoren auftreten (Auslöschung, $\ln$ bei 1, …).

Stabil umformen

Auslöschung lässt sich oft durch eine algebraisch äquivalente Umformung vermeiden. Standardtrick: mit der „dritten binomischen Formel“ erweitern.

$$D_1=\sqrt{a+h}-\sqrt a\;=\;\frac{(\sqrt{a+h}-\sqrt a)(\sqrt{a+h}+\sqrt a)}{\sqrt{a+h}+\sqrt a}\;=\;\frac{h}{\sqrt{a+h}+\sqrt a}=D_2.$$

Mathematisch gilt $D_1=D_2$. Numerisch subtrahiert $D_1$ zwei fast gleiche Zahlen (für kleines $h$). $D_2$ addiert nur positive Zahlen und dividiert, ist also stabil.

a =h = 10⁻ᵏ, gerechnet im Browser mit double (eps ≈ 1.1·10⁻¹⁶)

k	D1 = √(a+h) − √a	D2 = h / (√(a+h) + √a)	rel. Fehler von D1
1	4.880884817015163e-2	4.880884817015155e-2	1.7e-15
2	4.987562112088950e-3	4.987562112089027e-3	1.5e-14
3	4.998750624609638e-4	4.998750624609648e-4	2.0e-15
4	4.999875006239662e-5	4.999875006249610e-5	2.0e-12
5	4.999987500031722e-6	4.999987500062500e-6	6.2e-12
6	4.999998750587764e-7	4.999998750000624e-7	1.2e-10
7	4.999999880794803e-8	4.999999875000006e-8	1.2e-9
8	4.999999969612645e-9	4.999999987500000e-9	3.6e-9
9	5.000000413701855e-10	4.999999998750000e-10	8.3e-8
10	5.000000413701855e-11	4.999999999875000e-11	8.3e-8
11	5.000000413701855e-12	4.999999999987500e-12	8.3e-8
12	5.000444502911705e-13	4.999999999998750e-13	8.9e-5
13	4.996003610813204e-14	4.999999999999876e-14	8.0e-4
14	4.884981308350689e-15	4.999999999999988e-15	2.3e-2
15	4.440892098500626e-16	4.999999999999999e-16	1.1e-1
16	0	5.000000000000000e-17	1.0e+0
17	0	5.000000000000000e-18	1.0e+0

Beobachtung: $D_1$ verliert pro Zehnerpotenz von $h$ ungefähr eine gültige Stelle. Ab $h\lesssim\text{eps}\cdot a$ ist $\operatorname{rd}(a+h)=a$ und $D_1=0$. $D_2$ bleibt bis zum Schluss auf volle Genauigkeit.

Weitere typische Umformungen:

instabil	stabil	wann
$1-\cos x$	$2\sin^2\frac x2$	$x\approx0$
$\frac{-b+\sqrt{b^2-4ac}}{2a}$	$\frac{-2c}{b+\sqrt{b^2-4ac}}$	$b>0$, $b^2\gg4ac$
$e^x-1$	expm1(x)	$x\approx0$
$\ln(1+x)$	log1p(x)	$x\approx0$

(Die beiden unteren Zeilen sind Bibliotheksfunktionen, die es genau wegen dieser Auslöschung gibt. Ergänzung, nicht aus der Vorlesung.)

Übung 1

📄 Übungsblatt 1

Aufgabe 1: Archimedes' $\pi$ (Instabilität)

Kreis mit Radius $r=0.5$, also Umfang $\pi$. $s_n$ ist der Umfang des einbeschriebenen $n$-Ecks, $S_n$ der des umbeschriebenen. Rechnen in einfacher Genauigkeit (float, $\text{eps}\approx6\cdot10^{-8}$).

a) Herleitung. Eine Seite des $n$-Ecks ist $s_n/n=\sin(\pi/n)$ (Sehne bei Radius $\tfrac12$). Mit dem Halbwinkelsatz $\sin\frac\varphi2=\sqrt{\tfrac12(1-\cos\varphi)}$ und $\cos\varphi=\sqrt{1-\sin^2\varphi}$:

$$s_{2n}=2n\sin\frac{\pi}{2n}=2n\sqrt{\tfrac12\Big(1-\sqrt{1-(s_n/n)^2}\Big)}=n\sqrt{2\Big(1-\sqrt{1-(s_n/n)^2}\Big)},\qquad s_3=\tfrac32\sqrt3.$$

c) Warum instabil? Für großes $n$ ist $s_n/n$ klein, also $\sqrt{1-(s_n/n)^2}\approx1$. Die innere Klammer $1-\sqrt{\cdots}$ ist eine Auslöschung, und sie wird von Schritt zu Schritt schlimmer. Stabile Umformung wie oben mit $1-\sqrt{1-u^2}=\frac{u^2}{1+\sqrt{1-u^2}}$, $u=s_n/n$:

$$s_{2n}=n\sqrt{\frac{2u^2}{1+\sqrt{1-u^2}}}=s_n\sqrt{\frac{2}{1+\sqrt{1-(s_n/n)^2}}}.$$

Jetzt wird nur noch addiert, multipliziert und dividiert. Die Formel b) für $S_{2n}$ hat dieselbe Struktur (Plus im Nenner) und ist deshalb von Anfang an stabil.

Ergebnis in float (selbst nachgerechnet):

$n$	$s_n$ instabil	$s_n$ stabil
6	3.0000000	3.0000000
96	3.1410184	3.1410317
768	3.1430726	3.1415837
3072	3.1819804	3.1415915
12288	3.0000000	3.1415920
24576	4.2426405	3.1415920
49152	0.0000000	3.1415920

Die instabile Folge wird ab $n\approx200$ schlechter statt besser und bricht dann völlig zusammen ($1-\sqrt{1-u^2}$ wird in float exakt $0$).

Aufgabe 2: Integral-Rekursion (Fehlerfortpflanzung)

$$y_n=e^{-2}\int_0^1e^{2x}x^n\,dx,\qquad y_0=\tfrac12(1-e^{-2}),\qquad y_{n+1}=\tfrac12\big(1-(n+1)\,y_n\big).$$

a) Partielle Integration mit $u=x^{n+1}$, $v'=e^{2x}$: $\int_0^1e^{2x}x^{n+1}dx=\big[\tfrac12e^{2x}x^{n+1}\big]_0^1-\tfrac{n+1}2\int_0^1e^{2x}x^n\,dx$. Mit $e^{-2}$ multiplizieren ergibt die Rekursion.

b) Fehlerfortpflanzung. Ist $\tilde y_n=y_n+\Delta y_n$, dann

$\Delta y_{n+1}=-\frac{n+1}2\,\Delta y_n\quad\Rightarrow\quad|\Delta y_n|=\frac{n!}{2^n}\,|\Delta y_0|.$

Der Fehler wird pro Schritt mit $\frac{n+1}2$ multipliziert, ab $n\ge2$ also verstärkt, und zwar faktoriell. Schon die Rundung von $y_0$ ($|\Delta y_0|\approx10^{-17}$ in double) wächst bis $n=30$ auf $\frac{30!}{2^{30}}\cdot10^{-17}\approx10^{6}$. In double: $y_{20}\approx0.0435$ ist noch plausibel, $y_{25}\approx28.4$ und $y_{30}\approx-1.5\cdot10^7$ sind Unsinn (wahre Werte liegen in $(0,\tfrac12)$ und fallen).

c) Rückwärts ist stabil. Wegen $0\le y_n\le e^{-2}\int_0^1e^{2}x^n\,dx=\frac1{n+1}\to0$ ist $(y_n)$ eine Nullfolge. Umgestellt:

$y_n=\frac{1-2\,y_{n+1}}{n+1},\qquad\Delta y_n=-\frac{2}{n+1}\,\Delta y_{n+1}.$

Jetzt wird der Fehler pro Schritt durch $\frac{n+1}2$ geteilt. Der grobe Startfehler $y_{k+1}=0$ (statt $\approx\frac1{k+2}$) schrumpft bis $y_0$ um $\frac{2^{k+1}}{(k+1)!}$. Selbst nachgerechnet: $k=10$ gibt $y_0$ auf 5 Stellen, $k=20$ und $k=30$ auf volle double-Genauigkeit ($0.4323323583816936$).

Merke

Dieselbe Formel kann in eine Richtung instabil und in die andere stabil sein. Entscheidend ist, ob der Fehlerfaktor pro Schritt $>1$ oder $<1$ ist.