3 Fehler

3.1Fehlerarten¶

Bei numerischen Berechnungen spielen verschiedene Fehlerarten eine Rolle:

• Rundungsfehler (Zwischenergebnisse werden auf die darstellbare Genauigkeit gerundet).

• Datenfehler (die Eingabewerte entstehen aus (ungenauen) Messungen),

• Verfahrensfehler (die eingesetzten Algorithmen liefern auch bei exakter Rechnung ohne Rundungsfehler nur Näherungslösungen).

• Modellfehler (vereinfachte/idealisierte Darstellung der tatsächlichen physikalischen oder anderweitigen Gegebenheiten),

Dabei enthält der Datenfehler praktisch immer auch Rundungsfehler, da die Daten als Gleitkommazahlen an den Algorithmus übergeben werden müssen.

Obwohl die numerische Mathematik sich auf die (Vermeidung der) Auswirkungen von Rundungsfehlern konzentriert, spielen auch die anderen Fehlerarten und -quellen eine Rolle. Insbesondere zielt man gar nicht auf maximale Genauigkeit (was mit hohem Rechenaufwand verbunden ist), sondern nur darauf, dass die Auswirkungen von Verfahrens- und Rundungsfehlern nicht größer sind als zwangsläufige Abweichungen aus Daten- und Modellfehlern.

3.2Beispiel: Volumen der Erde¶

Wollen das Volumen der Erde aus gegebenem Radius $r$ (in Meter) mittels Kugelvolumen

berechnen.

Auftretende Fehler:

• Modellfehler: Die Erde ist keine exakte Kugel, sondern ein Ellipsoid (auch nicht exakt).

• Verfahrensfehler: Hier nicht vorhanden, sofern (bei exakten Daten) mit exaktem $\pi$ gerechnet wird.

• Datenfehler: Exakte Bestimmung von $r$ scheitert schon an der Struktur der Erdoberfläche (Berge, Täler). Selbst mit sehr präziser Messtechnik wird man nur ein Intervall $r\in[a,b]$ angeben können. Wählt man $r=\frac{a+b}{2}$, beträgt der Datenfehler $\frac{b-a}{2}$.

• Rundungsfehler: $r$ und $\pi$ müssen für die Darstellung mittels Gleitkommazahlen gerundet werden. Die Zahlen 4 und 3 sind exakt darstellbar. Die Auswirkungen der Multiplikationen und der Division auf die Genauigkeit des Ergebnisses werden weiter unten untersucht.

Aufgrund der unvermeidlichen großen Modell- und Datenfehler gibt es keinen Grund die Berechnungen mit sehr großer Genauigkeit auszuführen. Solange die Rundungsfehler und deren Auswirkungen auf die einzelnen Rechenoperationen kleiner als die Modell- und Datenfehler bleiben, ist die Genauigkeit groß genug. Hier kann sogar mit 3.14 statt $\pi$ gerechnet werden (Verfahrensfehler!), da die dadurch verursachten Abweichungen höchstens im Zentimeterbereich liegen.

3.3Bewertung von Algorithmen¶

Gegeben ist ein Algorithmus, der “ein Problem löst”. Genauer: Es gibt eine Menge $X\subseteq\bbR^d$ von möglichen Eingaben, eine Menge $Y\subseteq\bbR$ von möglichen Ausgaben und eine (exakte) Lösungsabbildung $f:X\to Y$.

Zentrale Frage: Wie hängt der Ausgabefehler vom Eingabefehler ab?

Dabei nehmen wir an, dass es eine unbekannte exakte Eingabe $x\in X$ gibt und eine bekannte mess- und rundungsfehlerbehaftete Eingabe $\tilde{x}$. Der relative Eingabefehler ist dann

Die numerische Auswertung von $f(x)$ erfolgt durch einen (potentiell fehlerhaften) Algorithmus $\tilde{f}:X\to Y$, der auf die fehlerbehafteten Daten $\tilde{x}$ angewendet wird. Für den relativen Ausgabefehler erhalten wir somit

Ziel einer Fehleranalyse sind Abschätzungen der Form

wobei

gelten soll. Die Zahl $\delta$ ist also eine obere Schranke an den zu erwartenden Datenfehler. Praktisch wird $\delta$ mindestens so groß sein wie das Maschinenepsilon.

Die Funktion $g$ ist meistens ein Monom $g(\delta)=c\,\delta^p$ mit $p>0$ und einer (hoffentlich kleinen) Konstante $c>0$. Je größer $p$, desto unempfindlicher reagiert der untersuchte Algorithmus auf Datenfehler. In jedem Fall muss $\lim_{\delta\to 0}g(\delta)=0$ gelten, damit die Fehlerabschätzung eine sinnvolle Aussage liefert (“Wenn wir den Datenfehler kleiner machen, wird auch der Lösungsfehler kleiner”).

Zentrale Werkzeuge für die Herleitung solcher Abschätzungen sind die Kondition der Abbildung $f$ und die Stabilität des Algorithmus $\tilde{f}$, welche im Folgenden eingeführt werden.

3.4Kondition¶

Einfaches Beispiel: Berechnung des Schnittpunktes zweier Geraden (IDVID 320).

• Verlaufen die Geraden nahezu senkrecht zueinander, dann gute Kondition.

• Verlaufen die Geraden nahezu parallel, dann schlechte Kondition.

Die Kondition hängt nicht nur von $f$ selbst, sondern auch von der konkreten Eingabe $x$ ab!

Die Kondition eines Problems bzw. der entsprechenden Abbildung $f$ kann man in Zahlen fassen. Für die konkrete Definition dieser Konditionszahl gibt es verschiedene Ansätze. Wir wählen zunächst den folgenden:

Der hier vorgestellte Konditionsbegriff führt immer zu linearen Fehlerzusammenhängen $g$ in (4). Dieser Ansatz ist in der Literatur auch als differentielle Fehleranalyse oder Störungstheorie 1. Ordnung zu finden.

Schlecht koniditionierte Probleme sind in der Praxis (leider) häufig anzutreffen, z.B. das bei der Computertomografie zu lösende mathematische Problem (siehe Veranstaltung zum wissenschaftlichen Rechnen).

3.5Beispiele zur Kondition¶

3.5.1Multiplikation¶

Betrachten

Für $x_1\neq 0$ und $x_2\neq 0$ erhalten wir $\kappa(x)=1$. Der Ausgabefehler ist als höchstens so groß wie der Eingabefehler.

Für $x_1=0$ oder $x_2=0$ können wir nur absolute Fehler betrachten:

Ist z.B. $x_1=0$ und $x_2$ ist sehr groß, also auch $\tilde{x}_2$ sehr groß, so kann das Produkt der gestörten Eingaben erheblich vom tatsächlichen Produkt 0 abweichen. Wähle z.B. $x_1=10^{-5}$ und $x_2=\tilde{x}_2=10^6$. Dann ist das berechnete Produkt 10 statt 0. Praktisch ist dieser Fall aber selten von Interesse, da die Null exakt als Gleitkommazahl darstellbar ist. Dennoch ist an dieser Stelle Vorsicht geboten!

a = 0.1 * 0.1 - 0.01
b = 10 ** 20
print(a * b)

173.4723475976807

3.5.2Division¶

Betrachten

Für $x_1\neq 0$ erhalten wir $\kappa(x)=1$. Der Ausgabefehler ist als höchstens so groß wie der Eingabefehler.

Für $x_1=0$ können wir nur absolute Fehler betrachten:

Die daraus zu ziehenden Schlüsse entsprechen denen bei der Multiplikation.

a = 0.1 * 0.1 - 0.01
b = 10 ** (-20)
print(a / b)

173.4723475976807

3.5.3Quadratwurzel¶

Betrachten

Für $x\neq 0$ erhalten wir $\kappa(x)=\frac{1}{2}$. Der Ausgabefehler ist also kleiner als der Eingabefehler.

Für $x=0$ erhalten wir den absoluten Fehler

der bei $\tilde{x}<1$ stets größer ist als der Eingabefehler $\tilde{x}$. Beachte hier, dass Eingabefehler $\tilde{x}<0$ zu undefiniertem Ergebnis bzw. komplexen Zahlen führen!

a = 0.01 - 0.1 * 0.1
print(a ** 0.5)

(8.064844237971058e-26+1.3170890159654386e-09j)

from numpy import sqrt

a = 0.01 - 0.1 * 0.1
print(sqrt(a))

nan

/tmp/ipykernel_25646/1240944823.py:4: RuntimeWarning: invalid value encountered in sqrt
  print(sqrt(a))

3.5.4Addition und Subtraktion¶

Betrachten

Für $x_1\neq 0$ und $x_2\neq 0$ erhalten wir

Haben $x_1$ und $x_2$ gleiches Vorzeichen, dann ist $\kappa(x)\leq 1$. Haben $x_1$ und $x_2$ jedoch unterschiedliches Vorzeichen, so kann $|x_1+x_2|$ klein sein, obwohl weder $|x_1|$ noch $|x_2|$ klein ist. In diesem Fall ist $\kappa$ sehr groß, die Addition als schlecht konditioniert! Diese Situation wird als Auslöschung bezeichnet.

Beispiel: Rechnen mit 5 Dezimalstellen, $x_1=0.120587$, $x_2=-0.120942$. Durch Runden der Eingabewerte entsteht ein Eingabefehler der Größenordnung 10^-4. Der Fehler im Ergebnis liegt jedoch bei 10^-2, ist also deutlich größer (Details: IDVID 340).

Weiteres Beispiel: Berechne die Differenz $a-b$, wobei statt $a$ und $b$ die auf 64-Bit-Gleitkommazahlen gerundeten Werte verwendet werden.

a = 1 + 11e-15
b = 1 + 10e-15
print(f'a = {a:.20f}')
print(f'b = {b:.20f}')

a = 1.00000000000001110223
b = 1.00000000000000999201

Mit exakten Werten erhalten wir $a-b=10^{-15}$. Der Eingabefehler ist kleiner als 10^-14, während der Ausgabefehler größer als 0.11 ist (IDVID 345). Somit weicht die tatsächliche Differenz um mehr als 11 Prozent vom berechneten Wert ab. Beachte hierbei: Wir haben exakt gerechnet, das Ergebnis also nicht gerundet. Ursache für den großen Fehler ist somit nicht ein konkreter Algorithmus (z.B. addiere exakt, runde dann auf 64-Bit-Gleitkommazahl), sondern das Problem selbst in Kombination mit leicht gestörten Eingabewerten. Berechnen wir die Differenz algorithmisch (also mit zusätzlichem Runden des Ergebnisses), so erhalten wir

c = a - b
print(f'a - b = {c}')
print( 'exakt = 1.0000000000000000e-15')

a - b = 1.1102230246251565e-15
exakt = 1.0000000000000000e-15

Nur die Ziffer vor dem Komma ist korrekt!

print(1 + c)         # Fehler in c nicht relevant
print(c * 10 ** 15)  # Fehler in c relevant

1.000000000000001
1.1102230246251565

3.5.5Sinus¶

Betrachten

Für $x\neq k\,\pi$, $k\in\bbZ$, erhalten wir

Für $x$ in der Nähe von Vielfachen von $\pi$ wird die Kondition also beliebig groß, während $f$ z.B. in der Nähe von $\frac{\pi}{2}$ sehr gut konditioniert ist.

from math import sin

def print_sin(a, b):

    sina = sin(a)
    sinb = sin(b)

    print(sina)
    print(sinb)

    print(f'Eingabefehler: {abs((a - b) / a)}')
    print(f'Ausgabefehler: {abs((sina - sinb) / sina)}')

print_sin(
    3.141593,
    3.141593123456789
)

-3.464102066193935e-07
-4.6986699585547026e-07
Eingabefehler: 3.929751219718377e-08
Ausgabefehler: 0.3563890060887287

Alle Mantissenziffern des für 3.141593 aus dem fehlerbehafteten Wert 3.141593123456789 berechneten Ergebnisses sind falsch, nur der Exponent stimmt. Der relative Fehler liegt bei 34 Prozent. Die Konditionszahl ist $\kappa(3.141593)\approx 9068997$. Beachte, dass die absoluten Fehler bei Eingabe und Ausgabe dennoch etwa gleich groß sind, das Ergebnis je nach Kontext also durchaus brauchbar sein kann.

print_sin(
    1.570796,
    1.570796123456789
)

0.9999999999999466
0.9999999999999793
Eingabefehler: 7.8595049270588e-08
Ausgabefehler: 3.2751579226443866e-14

Hier stimmen fast alle Mantissenziffern trotz gleichem absoluten Fehler in der Eingabe wie zuvor (der relative Eingabefehler ist hier sogar doppelt so groß). Der relative Ausgabefehler liegt im Bereich 10^-13, ist also im Größenordnungen kleiner als der relative Eingabefehler.

3.5.6Quadratische Gleichung¶

Betrachten das Lösen quadratischer Gleichungen $x^2+p\,x+q=0$, wobei wir nur die kleinere der beiden möglichen Lösungen berechnen möchten. Haben mit $D(f):=\{(p,q)\in\bbR^2:p^2>4\,q\}$ also

Für $x\neq 0$ und $f(p,q)\neq 0$ (also $p>0$ oder $q\neq 0$) erhalten wir

Für $q\leq 0$ gilt somit $\kappa(p,q)\leq 1$, also gute Kondition. Für $q>0$ wächst die Kondition mit $q$ und wird für $q\to\frac{p^2}{4}$ beliebig groß. Das Problem ist für große $q$ also sehr schlecht konditioniert (Details und grafische Interpretation: IDVID 350).

3.6Stabilität¶

Ein Algorithmus heißt also stabil, wenn er das tut, was man von ihm erwartet. Für die Formulierung präziserer Aussagen über den erwarteten Zusammenhang zwischen $\tilde{f}(x)$ und $f(x)$ gibt es verschiedene Ansätze. Wir gehen hier nur auf die beiden am weitesten verbreitetenen ein: Vorwärtsstabiliät und Rückwärtsstabilität.

Vorwärtsstabilität sichert also, dass der durch den verwendeten Algorithmus verursachte Fehler in der Größenordnung des auch bei exakter Rechnung nicht zu vermeidenden Fehlers liegt.

Für konkrete Algorithmen lassen sich Aussagen zur Vorwärtsstabilität oft nur schwer herleiten. Deshalb verwendet man fast immer den folgenden alternativen Stabilitätsbegriff.

Rückwärtsstabilität sichert also, dass die Ausgabe des Algorithmus der exakten Rechnung mit einer Eingabe entspricht, die im Rahmen der Größenordnung des ohnehin vorhandenen Eingabefehlers ebenfalls in Frage kommt.

Um Rückwärtsstabilität zu zeigen müssen zwei Fragen beantwortet werden:

1. Gibt es zu jeder Eingabe $x$ ein passendes $\hat{x}$?

2. Wie groß ist der Abstand zwischen $x$ und zugehörigem $\hat{x}$ höchstens?

Beweis: IDVID 360.

3.7Beispiele für Stabilitätsanalysen¶

3.7.1Addition/Subtraktion¶

Nach IEEE 754 erfolgt die Addition zweier Gleitkommazahlen nach folgendem Algorithmus:

Sind $x_1$, $x_2$ die beiden Summanden und ist $\varepsilon$ das Maschinenepsilon, so liefert der Algorithmus das Ergebnis

für ein $\varrho$ mit $|\varrho|\leq\varepsilon$.

Der einfache Algorithmus erlaubt eine Vorwärtsanalyse:

Die Kondition der Addition erfüllt $\kappa(x)\geq\frac{1}{2}$ (IDVID 370). Wir erhalten somit

Da der zu erwartende relative Eingabefehler praktisch immer größer als $2\,\varepsilon$ sein wird (Rundung jedes Summanden auf Gleitkommazahl), ist der untersuchte Algorithmus zur Addition als vorwärtsstabil anzusehen.

Zu Demonstrationszwecken untersuchen wir auch die Rückwärtsstabilität. Haben

mit $\hat{x}_1:=(1+\varrho)\,x_1$, $\hat{x}_2:=(1+\varrho)\,x_2$ also $\tilde{f}(x)=f(\hat{x})$. Der relative Fehler beim Ersetzen von $x$ durch $\hat{x}$ ist

Da die zu erwartenden Eingabefehler größer als $2\,\varepsilon$ sein werden, ist der Algorithmus als rückwärtsstabil anzusehen.

3.7.2Lange Summen¶

Die Summe von $n$ Zahlen soll berechnet werden:

Untersuchen dazu zunächst den einfachsten Algorithmus:

Mit Maschinenepsilon $\varepsilon$ und den aus den Rundungsschritten entstehenden Abweichungen $\varrho_2,\ldots,\varrho_n$ wird also

berechnet. Der Algorithmus ist für große $n$ nicht rückwärtsstabil (und auch nicht vorwärtstabil). Für die Konstante in der Definition der Rückwärtsstabilität erhalten wir $c=n$ (IDVID 375).

Betrachten nun folgenden alternativen Algorithmus für Summen mit $n=2^m$ Summanden:

In der letzten Iteration ist nur noch eine Summe mit zwei Summanden zu berechnen, welche das Endergebnis darstellt. Bezeichnen wir die bei den Rundungsschritten entstehenden Abweichungen mit $\varrho_{k,l}$, $k=0,\ldots,m-1$, $l=1,\ldots,2^{m-k-1}$, so entspricht der Algorithmus folgender Formel:

Für die Konstante in der Definition der Rückwärtsstabilität erhalten wir $c=m$ (IDVID 380). Der Algorithmus kann also auch für große $n$ noch als rückwärtsstabil (und somit auch vorwärtsstabil) angesehen werden, da $m=\log_2 n$ nur sehr langsam mit $n$ wächst.

3.7.3Quadratische Gleichung¶

Wir untersuchen zwei Algorithmen zur Berechnung der kleineren der beiden Lösungen der quadratischen Gleichung $x^2+p\,x+q=0$ auf Stabilität.

3.7.3.1Direktes Auswerten der Lösungsformel¶

Vorwärts- und auch Rückwärtsanalyse sind hier recht mühsam. Mindestens die Rückwärtsanalyse ist mit unseren Mitteln aber machbar. Wählen einen anderen Weg: Betrachten jeden Schritt separat und schließen aus er Kondition der Einzelschritte auf die (Nicht-)Stabilität des Gesamtalgorithmus.

• Schritt 1 ist stabil, liefert also bei fehlerfreier Eingabe $p$ einen Ausgabefehler unterhalb des Maschinenepsilons $\varepsilon$.

• Schritt 2 erhält eine fehlerbehaftete Eingabe, ist aber gut Konditioniert. Die Division ohne Runden wird also wieder einen Ausgabefehler in der Größenordnung von $\varepsilon$ liefern. Anschließendes Runden erhöht den relativen Fehler höchstens um den Faktor $1-\varepsilon$, sodass der Gesamtfehler nach Schritt 2 in der Größenordnung von $\varepsilon$ bleibt.

• Die Subtraktion in Schritt 3 ist schlecht konditioniert, wenn $q\approx\frac{p^2}{4}$, insbesondere für $q\leq 0$ aber gut konditioniert. Je nach Wert von $q$ kann der Ausgabefehler in Schritt 3 beliebig groß sein. Für $q\leq 0$ wird er bei $\varepsilon$ bleiben.

• Schritte 4 und 5 verhalten sich analog zu Schritt 2, also keine relevante Fehlerverstärkung.

• Schritt 6 ist schlecht konditioniert, wenn $p>0$ und $q\approx 0$; sonst gut konditioniert. Fehlerverhalten analog zu Schritt 3.

Der Algorithmus ist also nur stabil, wenn gilt:

3.7.3.2Umweg über Satz von Vieta¶

Wollen einen Algorithmus, der für $p\leq 0$ bei betragskleinem $q$ stabil ist. Der Satz von Vieta liefert $x_1\,x_2=q$ für die beiden Lösungen der quadratischen Gleichung. Wir können also zunächst die größere Lösung berechnen und dann daraus die kleinere. Die ersten 5 Schritte bleiben unverändert:

Für $p\leq 0$ ist Schritt 6 nun gut Konditioniert (unabhängig vom Wert von $q$). Auch Schritt 7 ist gut konditioniert solange das Ergebnis aus Schritt 6 nicht zu nah bei 0 liegt. Allerdings ist Schritt 6 für $p>0$ und $q\approx 0$ nun schlecht konditioniert. Der Algorithmus ist also nur stabil, wenn gilt:

3.7.3.3Kombinierter Algorithmus¶

Um die kleinere Lösung einer quadratischen Gleichung stabil zu erhalten, oder auch um beide Lösungen zu erhalten, sollten also erst die obigen Schritte 1 bis 5 ausgeführt werden. Anschließend ist zu unterscheiden:

• Wenn $p>0$, dann Subtrahiere die Wurzel von $-\frac{p}{2}$.

• Wenn $p<0$, dann Addiere die Wurzel zu $-\frac{p}{2}$.

Anschließend kann über die Beziehung $x_1\,x_2=q$ (Vieta!) die andere Lösung berechnet werden.

Dieser Algorithmus ist stabil, wenn $q\ll\frac{p^2}{4}$.

Die bei $q\approx\frac{p^2}{4}$ auftretende Instabilität aufgrund von Auslöschungseffekten spielt praktisch kaum eine Rolle. In diesem Fall sind beide Lösungen sehr nah beieinander. Sind beide Lösungen hinreichend weit von der Null entfernt, werden diese trotzdem korrekt berechnet, da dann nur der durch Auslöschung entstandene absolute Fehler relevant ist für den relativen Fehler in den Lösungen. Und dieser ist sehr klein. Liegen jedoch beide Lösungen nah bei Null, so überträgt sich der durch Auslöschung entstandene große relative Fehler auf den relativen Fehler der Lösungen.

3.8Gesamtfehler¶

Faustregel: Gute Kondition des Problems $f$ und Stabilität des Algorithmus $\tilde{f}$ garantieren einen hinnehmbar kleinen Gesamtfehler. Schelchte Kondition oder fehlende Stabilität lassen keine brauchbaren Ergebnisse erwarten.

Beweis: IDVID 390.

Die Abschätzung zeigt, dass der Ausgabefehler um so kleiner wird, je kleiner der Eingabefehler ist (wenn $\kappa(\tilde{x})$ bei $\tilde{x}\to x$ beschränkt bleibt).