Aus Wikiversity

Inhaltsverzeichnis

1 1.1 b-adische Brüche
2 1.2 Rechnen auf einem Computer
3 1.3 Stabilität und Kondition
4 1.4 Differentielle Fehleranalyse
- 4.1 Beispiel 1.11

[Bearbeiten] 1.1 b-adische Brüche

Alle Berechnungen auf digitalen Rechenanlagen können grundsätzlich nur mit endlicher Genauigkeit durchgeführt werden. Deshalb spielt die verwendete Zahlendarstellung eine wichtige Rolle.

Bekanntlich ist die Zahl $x := 311.57 \ldots$ im Dezimalsystem gleichbedeutend mit

$3 \cdot 10^2 + 1 \cdot 10^1 + 1 \cdot 10^0 + 5 \cdot 10^{-1} + 7 \cdot 10^{-2} + \ldots$

Man spricht in diesem Fall auch von einem $b$ -adischen Bruch zur Basis $b = 10$ . Statt der bekannten Basis $b : = 10$ in unserem Dezimalsystem kann man auch eine andere Basis $b \in \mathbb N$ mit $b \ge 2$ wählen. Wäre $b : = 8$ , so hätte beispielsweise die Zahl mit der Ziffernfolge $311.57 \ldots$ im Dezimalsystem den Wert

$3 \cdot 8^2 + 1 \cdot 8^1 + 1 \cdot 8^0 + 5 \cdot 8^{-1} + 7 \cdot 8^{-2} + \ldots$

Sind $m_i \in \{0, 1, \ldots, b - 1\}$ insbesondere die Stellen einer Zahl bezüglich der Basis $b$ , so entspricht diese also im Dezimalsystem der reellen Zahl

$x = \pm \Bigl( m_nb^n + m_{n-1}b^{n-1} + \ldots + m_0b^0 + m_{-1}b^{-1} + m_{-2}b^{-2} + \ldots \Bigr) = \pm \sum^n_{i=-\infty} m_ib^i.$

Praktisch wichtig sind dabei die Basen

$b = 2: \quad \text{Dualbasis}, \quad m_i \in \{0, 1\},$

$b = 10: \quad \text{Dezimalbasis}, \quad m_i \in \{0, 1, \ldots, 9\},$

$b = 16: \quad \text{Hexadezimalbasis}, \quad m_i \in \{0, 1, \ldots, 9, A, B, C, D, E, F\}.$

Die Babylonier beispielsweise verwendeten die Basis $b = 60$ (Sexagesimalsystem).

[Bearbeiten] Beispiel 1.1

(1) $b = 2.$ Dem sog. dyadischen Bruch $11100.10$ entspricht im Dezimalsystem die Zahl

$1 \cdot 2^4 + 1 \cdot 2^3 + 1 \cdot 2^2 + 0 \cdot 2^1 + 0 \cdot 2^0 + 1 \cdot 2^{-1} = 16 + 8 + 4 + \frac{1}{2} = 28.5.$

(2) $b = 16.$ Der $b$ -adische Bruch $C 17. E$ zur Basis $b = 16$ bedeutet im Dezimalsystem die Zahl

$12 \cdot 16^2 + 1 \cdot 16^1 + 7 \cdot 160 + 14 \cdot 16^{-1} = 3095.875.$

Will man umgekehrt eine Dezimalzahl bezüglich einer anderen Basis als $b : = 10$ beschreiben, so verfährt man wie folgt.

[Bearbeiten] Beispiel 1.2

Gegeben sei die Zahl $x : = 60.125$ im Dezimalsystem. Diese sei nun mittels einer anderen Basis dargestellt.
(1) $b = 2.$ (Als Vielfache der Potenzen von 2 stehen nur die Ziffern 0 und 1, für die erste Stelle nur die 1 zur Verfügung.) Man ermittelt zuerst die höchste Potenz von 2, die in 60 steckt. Dies ist wegen $25 = 32$ und $26 = 64$ die 5. Weiter enthält $60.125 - 2 5 = 28.125$ das 1-fache von $24 = 16$ , die Zahl $28.125 - 16 = 12.125$ das 1-fache von $23 = 8$ , die Zahl $4.125$ das 1-fache von $22 = 4$ , die Zahl $0.125$ das 0-fache von $21 = 2$ sowie das 0-fache von $20 = 1$ , die Zahl $0.125$ das 0-fache von $2 - 1 = 0.5$ sowie das 0-fache von $2 - 2 = 0.25$ und schließlich das 1-fache von $2 - 3 = 0.125$ , so dass man zu folgender Darstellung gelangt:

111100.001.

(2) $b = 8.$ (Als Vielfache der Potenzen von 8 stehen jetzt die Ziffern $0, 1, \ldots, 7$ zur Verfügung bzw. für die erste Stelle die Ziffern $1, 2, \ldots, 7$ .) Es ergibt sich

$60.125 = 7 \cdot 8^1 + 4 \cdot 8^0 + 1 \cdot 8^{-1}.$

Also entspricht $60.125$ im Dezimalsystem der Zahl $74.1$ zur Basis $b = 8$ .

(3) $b = 16.$ Man erhält

$60.125 = 3 \cdot 16^1 + 12 \cdot 16^0 + 2 \cdot 16^{-1}.$

Es ergibt sich zur Basis $b = 16$ die Zahl $3 C .2$ .

[Bearbeiten] Satz 1.3

Sei $b \in \mathbb N$ und $b \ge 2$ . Dann lässt sich jede reelle Zahl in einen $b$ -adischen Bruch zur Basis $b$ entwickeln.

[Bearbeiten] 1.2 Rechnen auf einem Computer

Wir gehen nun von einer Zahlendarstellung von $x \in \mathbb R$ mittels der Basis $b \in \mathbb N$ mit $b \ge 2$ und den Ziffern $m_i \in \{0, 1, \ldots, b - 1\}$ aus, d. h. von einer Darstellung

$x = \pm \sum^n_{i=-\infty} m_ib^i$

mit $m_n \neq 0$ . (Für $b > 10$ kann $m i$ also auch eine Ziffer mit mehr als einer Stelle sein.) Durch Abschneiden dieses unendlichen Ausdrucks ergibt sich eine endliche Zahlendarstellung

$x = \pm \Bigl( m_nb^n + m_{n-1}b^{n-1} + \ldots + m_0b^0 + m_{-1}b^{-1} + \ldots + m_{-p}b^{-p} \Bigr) = \pm \sum^n_{i=-p} m_ib^i.$

Digitale Rechenanlagen, kurz Computer oder Rechner, arbeiten meist mit einer normalisierten (endlichen) Gleitkomma-Darstellung reeller Zahlen

(1.1) $a := \pm M \cdot b^E,$

wobei die sog. Mantisse $M$ von der Form

$M = 0.\alpha_1 \ldots \alpha_r$

und der Exponent $E$ von der Form

$E = \pm e_1 \ldots e_s$

mit Ziffern $\alpha_i, e_j \in \{0, \ldots, b - 1\}$ sind und insbesondere $\alpha_1 \neq 0$ gilt. Diese Darstellung ist wegen $\alpha_1 \neq 0$ eindeutig. Für die Darstellung der Null setzt man $M : = 0$ und $E$ beliebig. Die Speicherung von $a$ in (1.1) erfolgt in der Form

$a \to (\pm) [\alpha_1 \ldots \alpha_r] (\pm) [e_1 \ldots e_s]$

und benötigt

r

Ziffern + 1 Vorzeichen für die Mantisse,

s

Ziffern + 1 Vorzeichen für den Exponenten.

[Bearbeiten] Beispiel 1.4

Sei $b : = 10$ .
(1) Die Zahl $- 30.421$ lautet (bei Nichtberücksichtigung der Größen $r$ und $s$ ) in normalisierter Gleitkomma-Darstellung $- 0.30421 \cdot 10^2$ . Letztere Darstellung schreibt man z.B. für $r : = 6$ und $s : = 2$ auch in der Form $- 0.304210 E + 02$ oder $- 0.304210 10 + 02$ .

(2) Die Zahl $- 0.00030421$ lautet in der normalisierten Gleitkomma-Darstellung für $r : = 6$ und $s : = 2$ z. B. $- 0.304210 E - 03$ oder $- 0.304210 10 - 03$ .

Eine normalisierte Gleitkomma-Darstellung mit der Basis $b$ , beispielsweise $b : = 10$ oder $b : = 2$ , bestimmt die Menge $A : = A (b; r; s)$ reeller Zahlen, die auf dem Rechner exakt dargestellt werden können, die sog. Maschinenzahlen. Eine solche Zahlendarstellung ermöglicht also nur die Repräsentation einer endlichen Teilmenge der reellen Zahlen. Und zwar ist offenbar insbesondere die kleinste darstellbare positive Zahl $a min$ durch

$M := 0.1, \quad E := -\underbrace{(b - 1)(b - 1) \ldots (b - 1)}_{s-mal}$

und die größte positive Zahl $a max$ durch

$M := 0.\underbrace{(b - 1)(b - 1) \ldots (b - 1)}_{r-mal}, \quad E := +\underbrace{(b - 1)(b - 1) \ldots (b - 1)}_{s-mal}$

gegeben. Die Mantisse $M$ von $a min$ entspricht offenbar der Dezimalzahl $b - 1$ und die von $a max$ der Dezimalzahl

$(b - 1) \sum^r_{i=1} b^{-i} = (b - 1) \left[ \frac{1 - b^{-r-1}}{1 - b^{-1}} - 1 \right] = \frac{b (b^{r+1} - 1)}{b^{r+1}} - b + 1 = 1 - b^{-r}.$

Der Exponent $E$ für beide Zahlen ist bis auf das Vorzeichen gegeben durch die Dezimalzahl

$(b - 1) \sum^{s-1}_{i=0} b^i = (b - 1) \frac{1 - b^s}{1 - b} = b^s - 1.$

Somit haben $a min$ und $a max$ den Dezimalwert

$a_\min = b^{-b^s}, \quad a_\max = (1 - b^{-r})b^{b^s-1}.$

Ist $D$ die reelle Zahlenmenge

$D := [-a_\max, - a_\min] \cup \{0\} \cup [a_\min, a_\max],$

so können also insbesondere Zahlen $x \notin D$ nicht auf dem Rechner wiedergegeben werden. Im Fall $x > a max$ und $x < - a max$ melden alle Rechner normalerweise einen Exponentenüberlauf („overflow“), während sie im Fall $x \in (-a_\min, a_\min)$ meist keine Meldung machen und $x = 0$ setzen.

Ferner ist offenbar nicht jede Zahl $x \in D$ auf dem Rechner, d. h. als $x \in A$ darstellbar (z. B. trifft dies für die Zahlen $π$ und $\sqrt{2}$ zu). Somit stellt sich das Problem, eine Zahl $x \in D$ durch eine Zahl aus $A$ zu approximieren. Man verwendet hierzu einen Rundungsoperator $\operatorname{rd}: D \to A$ , der jeder Zahl $x \in D$ eine Zahl $\operatorname{rd}(x) \in A$ zuordnet, welche sinnvollerweise der folgenden Beziehung genügt:

(1.2) $|x - \operatorname{rd}(x)| = \min_{a \in A} |x - a|.$

Im Fall, dass die Aufgabe in (1.2) zwei Lösungen besitzt, rundet man dabei normalerweise (wir legen dies hier auch so fest) z. B. für $b : = 10$ und eine Endziffer 5 nach oben.

[Bearbeiten] Beispiel 1.5

Sei $b : = 10, r : = 4$ und $s : = 2$ . Dann gilt

$\operatorname{rd}(3.14159) = 0.3142_{10} + 01,$

$\operatorname{rd}(14.2842) = 0.1428_{10} + 02,$

$\operatorname{rd}(0.142749) = 0.1427_{10} + 00,$

$\operatorname{rd}(0.14275) = 0.1428_{10} + 00.$

Allgemein kann man für $b : = 10$ und eine Mantissenlänge $r$ die zu einer beliebigen Zahl $x \in D$ gehörende Maschinenzahl $\operatorname{rd}(x) \in A$ folgendermaßen finden. Es sei dazu $x \in D$ zunächst in der Form $x = a \cdot 10^q$ dargestellt, wobei $q \in \mathbb Z$ und

$|a| = 0.\alpha_1 \alpha_2 \ldots \alpha_r \alpha_{r+1} \ldots$

mit $0 \le \alpha_i \le 9$ und $\alpha_1 \neq 0$ seien. Insbesondere ist also $|a| \ge 0.1$ . Zu $a$ bildet man nun

$\tilde a :=\begin{cases} 0.\alpha_1 \alpha_2 \ldots \alpha_r, & \text{falls } 0 \le \alpha_{r+1} \le 4, \\ 0.\alpha_1 \alpha_2 \ldots \alpha_r + 10^{- r}, & \text{falls } \alpha_{r+1} \ge 5 \end{cases}$

und setzt dann

$\operatorname{rd}(x) := \sgn(x) \cdot \tilde a \cdot 10^q.$

Offenbar ist die Zahl $\operatorname{rd}(x)$ für jedes $x \in D$ eine Maschinenzahl, d. h. $\operatorname{rd}(x) \in A$ .

[Bearbeiten] Beispiel 1.6

Für $b : = 10, r : = 4$ und $s : = 2$ folgt

$\operatorname{rd}(0.99997_{10} + 98) = 0.1000_{10} + 99,$

$\operatorname{rd}(0.012345_{10} - 9) = 0.1235_{10} - 10.$

Für den mit der Rundung verbundenen absoluten Fehler hat man

$|x - \operatorname{rd}(x)| \le 5 \cdot 10^{-(r+1)} \cdot 10^q = \frac{1}{2} 10^{-r} 10^q$

und für den relativen Fehler, sofern $x \neq 0$ ist,

$\frac{|x - \operatorname{rd}(x)|}{|x|} \le \frac{0.5 \cdot 10^{-r} 10^q}{|a| 10^q} \le \frac{0.5 \cdot 10^{-r}}{0.1} = \frac{1}{2} 10^{-r+1}.$

Bei einer analogen Definition der Rundungsoperation für die Basis $b$ erhält man

$|x - \operatorname{rd}(x)| \le \frac{1}{2} b^{- r}b^q, \quad \frac{|x - \operatorname{rd}(x)|}{|x|} \le \frac{1}{2} b^{-r+1}.$

Diese Vorgehensweise führt auf die Definition der sogenannten relativen Maschinengenauigkeit

$eps := \frac{1}{2} b^{-r+1}.$

Mit dieser gilt also

(1.3) $\operatorname{rd}(x) = x (1 + \varepsilon), \quad |\varepsilon| \le eps,$

wie man mit der Setzung $\varepsilon := -(x - \operatorname{rd}(x))/x$ sieht.

[Bearbeiten] Beispiel 1.7 (IEEE-Standard)

$\begin{array}{l|c|c} & \text{single precision} & \text{double precision} \\ \hline a_\min & 1.10 \cdot 10^{- 38} & 2.23 \cdot 10^{- 308} \\ a_\max & 3.40 \cdot 10^{+ 38} & 1.80 \cdot 10^{+ 308} \\ eps & 0.60 \cdot 10^{- 7} & 1.11 \cdot 10^{- 16} \end{array}$

Die arithmetischen Grundoperationen $+ , - , * , /$ werden auf digitalen Rechnern durch sog. Gleitpunktoperationen ersetzt, welche Maschinenzahlen wieder auf Maschinenzahlen abbilden. Sind $a, b \in A$ , ist „ $\circ$ “ eine der vier Grundoperationen, $c := a \circ b$ und $c \in D$ , so definiert man die zugehörige Gleitpunktoperation $gl(a \circ b)$ durch

$gl(a \circ b) := \operatorname{rd}(a \circ b).$

Für sie gilt nach (1.3)

$gl(a \circ b) = (a \circ b) (1 + \varepsilon), \quad |\varepsilon| \le eps.$

Für das Ergebnis $c := a \circ b$ kann natürlich auch $c \notin D$ gelten. In diesem Fall meldet der Rechner normalerweise einen Exponentenüberlauf oder setzt er $c : = 0$ .

[Bearbeiten] 1.3 Stabilität und Kondition

Unter einem Algorithmus verstehen wir eine eindeutig festgelegte Folge von elementaren Operationen $+ , - , * , /$ , d. h. ein eindeutig festgelegtes Verfahren zur numerischen Lösung eines Problems oder einer Klasse von Problemen. Zumeist gibt es mehrere unterschiedliche Algorithmen zur Lösung eines Problems, die bei exakter Rechnung zu demselben Ergebnis führen würden, die dies jedoch numerisch aufgrund von Rundungsfehlern, Speicherplatz etc. auf dem Computer im Allgemeinen nicht tun.

[Bearbeiten] Beispiel 1.8

Seien $a = 0.03615, b = 62.88$ und $c = 62.73$ . Dann gilt bei exakter Rechnung

a + b - c = 0.18615.

Bei 4-stelliger Rechnung ergibt sich hingegen

g l (a + g l (b - c)) = g l (0.03615 + 0.1500) = 0.1862

mit 4 in Bezug auf das exakte Ergebnis korrekten Stellen oder alternativ

g l (g l (a + b) - c) = g l (62.92 - 62.73) = 0.1900

mit einem Ergebnis, welches lediglich 2 korrekte Stellen aufweist.

Die $g l$ -Operationen genügen also weder dem Assoziativ- noch dem Distributivgesetz!

Das Ziel der numerischen Mathematik ist die Konstruktion und Analyse von effizienten Algorithmen. Dabei meinen wir mit effizient, dass sie

(a) schnell und sparsam sein, d. h. möglichst wenige Rechenoperationen und damit bei geringem Speicherplatzbedarf wenig Rechenzeit zur Berechnung des gewünschten Resultates benötigen sollten,

(b) numerisch stabil sein sollten, d. h. die Verfälschung der Resultate durch Rundungsfehler nicht wesentlich größer sein sollte als die Verfälschung durch die Eingabefehler.

Unter Eingabefehler versteht man dabei die durch Rundung der Eingabedaten auf dem Rechner sich bereits ergebenden (im Allgemeinen kleinen) Fehler. Ein mathematisches Problem kann aber selbst schon „gutartig“ oder „bösartig“ sein. So nennt man ein mathematisches Problem gut konditioniert, wenn kleine Störungen in den Daten auch nur kleine Fehler in den Lösungen zur Folge haben, und es heißt schlecht konditioniert (ill-conditioned), wenn kleine Störungen in den Daten große Fehler in den Lösungen bewirken können. Bei schlecht konditionierten Problemen ziehen also Eingabefehler auf dem Rechner üblicherweise automatisch große Fehler in den erzielten Resultaten nach sich und zwar für jeden Algorithmus.

[Bearbeiten] Beispiel 1.9

(1) Die Subtraktion $c : = a - b$ etwa gleich großer reeller Zahlen $a$ und $b$ ist ein schlecht konditioniertes Problem. Denn sind $a$ und $b$ mit Fehlern $\varepsilon_a$ und $\varepsilon_b$ gestört, d. h. hat man statt $a$ und $b$

$\tilde a = a (1 + \varepsilon_a), \quad \tilde b = b (1 + \varepsilon_b),$

dann ergibt sich

$\tilde c = \tilde a - \tilde b = a (1 + \varepsilon_a) - b (1 + \varepsilon_b) = (a - b) \left( 1 + \frac{a}{a - b} \varepsilon_a - \frac{b}{a - b} \varepsilon_b \right).$

Für $a \approx b$ (und $|a - b| \ll |a|$ bzw. $|a - b| \ll |b|$ ) hat man bezüglich $c : = a - b$ im Ergebnis eine Fehlerverstärkung gegenüber $\varepsilon_a$ bzw. $\varepsilon_b$ von

$\left| \frac{a}{a - b} \right| \gg 1, \quad \left| \frac{b}{a - b} \right| \gg 1$

Auf dem Rechner führt dies zum Phänomen der Auslöschung von korrekten Stellen. Hat man z. B. bei 8-stelliger Rechnung

$\tilde a = 0.123\, 467\, xx$ (6 korrekte Stellen),

$\tilde b = 0.123\, 456\, 1x$ (7 korrekte Stellen),

so erhält man

$\tilde a - \tilde b = 0.000\, 011\, xx = 0.11x\, x00\, 00_{10} - 4$ (2 korrekte Stellen).

(2) Das Problem, den Schnittpunkt zweier Geraden $g$ und $h$ , die fast parallel zueinander sind, zu berechnen, ist schlecht konditioniert. Um dies einzusehen, mache man sich geometrisch klar, dass im Fall zweier sich im rechten Winkel schneidender Geraden kleine „Störungen“ der Geraden auch nur kleine Störungen bezüglich des gemeinsamen Schnittpunktes zur Folge haben, während solche Störungen großen Einfluss haben, wenn die Geraden fast parallel zueinander sind.

Bei der Konstruktion von Algorithmen sollte man also, wenn möglich, schlecht konditionierte Teilprobleme vermeiden.

[Bearbeiten] Beispiel 1.10

(1) Der Ausdruck

c : = a 2 - b 2

sollte mittels der numerisch stabileren Darstellung

c : = (a - b)(a + b)

berechnet werden.

(2) Die Funktion

$f(x) := \frac{1}{1 + 2x} - \frac{1 - x}{1 + x}$

erlaubt die stabilisierte Darstellung

$y = \frac{2x^2}{(1 + x)(1 + 2x)},$

welche sich zur Berechnung von Funktionswerten für $|x| \ll 1$ anbietet. Man mache sich klar, dass bei der Auswertung der stabilisierten Darstellung keine Auslöschung auftritt.

[Bearbeiten] 1.4 Differentielle Fehleranalyse

Der Einfluss von Störungen in den Daten auf die Lösung eines Problems sowie die Fortpflanzung von Eingangs- und Rundungsfehlern bei numerischen Algorithmen kann durch die sogenannte differentielle Fehleranalyse untersucht werden. Zu deren Beschreibung nehmen wir an, dass ein Problem bzw. eine Berechnungsvorschrift mittels einer zweimal stetig differenzierbaren Funktion $f: \mathbb R^n \to \mathbb R^m$ durch die Gleichung

f (x) = y

bzw. gleichbedeutend durch die Gleichungen

(1.4) $f_i(x) = y_i, \quad i = 1, \ldots, m$

beschrieben wird. Dabei ist also $x \in \mathbb R^n$ der Daten- und $y \in \mathbb R^m$ der Ergebnis-Vektor. Für

f (x + Δ x) = : y + Δ y

gilt dann nach dem Satz von Taylor zeilenweise

(1.5) $\Delta y_i = f_i(x + \Delta x) - f_i(x) = \sum^n_{j=1} \frac{\partial f_i}{\partial x_j} (x) \Delta x_j + \mathcal O \left( \max_{j=1, \ldots, n} |\Delta x_j|^2 \right), \quad i = 1, \ldots, m,$

so dass für hinreichend kleine $| Δ x j |$ der Restterm $\mathcal O \left( \max_{j=1, \ldots, n} |\Delta x_j|^2 \right)$ vernachlässigt werden kann und folglich der dominierende relative Fehler gegeben ist durch

$\frac{|\Delta y_i|}{|y_i|} \approx \sum^n_{j=1} \left| \frac{\partial f_i}{\partial x_j} (x) \frac{x_j}{f_i(x)} \right| \left| \frac{\Delta x_j}{x_j} \right| = \sum^n_{j=1} k_{ij}(x) \left| \frac{\Delta x_j}{x_j} \right|, \quad i = 1, \ldots, m$

mit

$k_{ij}(x) :=\left| \frac{\partial f_i}{\partial x_j} (x) \right| \left| \frac{x_j}{f_i(x)} \right|$

Die Größen $k i j (x)$ entscheiden demnach über den Einfluss der relativen Fehler $| Δ x j | / | x j |$ in den Daten auf den relativen Fehler $| Δ y i | / | y i |$ im Ergebnis. Sie werden deshalb häufig auch Verstärkungsfaktoren genannt. Im Fall, dass die Ausgangsgleichung einen Algorithmus beschreibt, sagt man, dass dieser stabil ist, wenn alle $k i j (x)$ „klein“, idealerweise ungefähr gleich 1 sind. Anderenfalls sagt man, er ist instabil.

Im Fall, dass die Ausgangsgleichung ein mathematisches Problem beschreibt, spricht man bei den $k i j (x)$ auch von Konditionszahlen (engl. to condition = bedingen, bestimmen). Sind die Konditionszahlen dem Betrag nach groß, hat man also ein schlecht konditioniertes, anderenfalls ein gut konditioniertes Problem. Für manche Zwecke ist aber diese Definition von Konditionszahlen unpraktisch, so dass auch andere Größen als Konditionszahlen bezeichnet werden (vgl. Definition 2.18).

[Bearbeiten] Beispiel 1.11

Die Lösungen $λ 1$ und $λ 2$ einer quadratischen Gleichung

(1.6)

x 2 - 2 p x + q = 0

sind gegeben durch

(1.7) $\lambda_{1, 2} := p \pm \sqrt{p^2 - q},$

wobei wir hier davon ausgehen, dass die Gleichung zwei unterschiedliche reelle Lösungen besitzt, also

p 2 - q > 0

ist. Zur Analyse der Fehlerempfindlichkeit der beiden Lösungen von (1.6) in Abhängigkeit von den Eingabedaten $p$ und $q$ betrachten wir diese nun als Funktionen von $p$ und $q$ . Wir untersuchen dazu die beiden Gleichungen

(1.8) $\lambda_1(p, q) = p + \sqrt{p^2 - q}, \quad \lambda_2(p, q) = p - \sqrt{p^2 - q}.$

(Im Vergleich mit (1.4) ist $x : = (p, q), m = 2, f i : = λ i$ und entsprechen die rechten Seiten in (1.8) den $y i$ .) Man hat dafür

(1.9) $\lambda_1 + \lambda_2 = 2p, \quad \lambda_1 - \lambda_2 = 2\sqrt{p^2 - q}, \quad \lambda_1 \lambda_2 = q.$

Damit errechnet man

$\frac{\partial \lambda_{1, 2}}{\partial p} = 1 \pm \frac{p}{\sqrt{p^2 - q}} = \frac{\sqrt{p^2 - q} \pm p}{\sqrt{p^2 - q}} = \pm \frac{\lambda_{1, 2}}{\lambda_1 - \lambda_2},$

$\frac{\partial \lambda_{1, 2}}{\partial p} = \mp \frac{1}{2\sqrt{p^2 - q}} = \mp \frac{1}{\lambda_1 - \lambda_2}.$

Hieraus ergeben sich die Verstärkungsfaktoren

$k_{11} := \left| \frac{\partial \lambda_1}{\partial p} \frac{p}{\lambda_1} \right| = \left| \left( \frac{2\lambda_1}{\lambda_1 - \lambda_2} \right) \left( \frac{\lambda_1 + \lambda_2}{2\lambda_1} \right) \right| = \frac{|1 + (\lambda_1/\lambda_2)|}{|1 - (\lambda_1/\lambda_2)|},$

$k_{12} := \left| \frac{\partial \lambda_1}{\partial q} \frac{q}{\lambda_1} \right| = \left| \left( \frac{1}{\lambda_1 - \lambda_2} \right) \left( \frac{\lambda_1 \lambda_2}{\lambda_1} \right) \right| = \frac{1}{|1 - (\lambda_1/\lambda_2)|},$

$k_{21} := \left| \frac{\partial \lambda_2}{\partial p} \frac{p}{\lambda_2} \right| = \left| \left( \frac{2\lambda_2}{\lambda_1 - \lambda_2} \right) \left( \frac{\lambda_1 + \lambda_2}{2\lambda_2} \right) \right| = k_{11},$

$k_{22} := \left| \frac{\partial \lambda_2}{\partial q} \frac{q}{\lambda_2} \right| = \left| \left( \frac{1}{\lambda_1 - \lambda_2} \right) \left( \frac{\lambda_1 \lambda_2}{\lambda_2} \right) \right| = k_{12}.$

Die Bestimmung der Lösungen von (1.6) ist somit für $\lambda_1 \approx \lambda_2$ ein schlecht konditioniertes Problem. Zur Veranschaulichung geben wir ein Zahlenbeispiel. Es seien $p : = 2$ und $q : = 3.999$ . Dann sind $λ 1 = 2.01$ und $λ 2 = 1.99$ die beiden Nullstellen von (1.6). Für die Verstärkungsfaktoren ergibt sich in diesem Fall

$k_{11} = k_{21} = \frac{|1 + (\lambda_1/\lambda_2)|}{|1 - (\lambda_1/\lambda_2)|} \approx 200, \quad k_{22} = k_{12} = \frac{1}{|1 - (\lambda_1/\lambda_2)|} \approx 99.5.$

Somit ist zu erwarten, dass Eingabefehler in den Daten $p$ und $q$ in Bezug auf die Lösungen $λ 1$ und $λ 2$ von (1.6) um den 100- bis 200-fachen Wert verstärkt werden.

Wir wollen nun zwei unterschiedliche Algorithmen zur Berechnung von

$\lambda_1 := p + \sqrt{p^2 - q}, \quad \lambda_2 := p - \sqrt{p^2 - q}$

betrachten und zwar unter den Bedingungen

(1.10) $p^2 - q > 0, \quad |q| \ll p^2, \quad p < 0.$

In diesem Fall ist offenbar

$\lambda_1 \approx -|p| + |p| = 0, \quad \lambda_2 \approx - |p| - |p| = -2 |p|,$

d. h. für nicht zu kleine $| p |$ auch $\lambda_1 \gg \lambda_2$ und somit das Problem der Bestimmung der Lösungen der quadratischen Gleichung (1.6) gut konditioniert. Ein „Lösungsalgorithmus“ könnte nun zunächst darin bestehen, hintereinander die folgenden Größen zu berechnen

$u := p^2, \quad v := u - q, \quad w := \sqrt{v} \ge 0.$

Wegen $p < 0$ sollte man als nächstes den unkritischen Wert

λ 2 : = p - w

berechnen. Zur Berechnung von $λ 1$ betrachten wir nun zwei Varianten (vgl. (1.9)):

$\text{Variante A}: \quad \lambda_1 := p + w,$

$\text{Variante B}: \quad \lambda_1 := q/\lambda_2.$

Da unter den Voraussetzungen (1.10) $w \approx - p$ gilt, tritt bei Variante A zwangsläufig Auslöschung auf. Betrachtet man $λ 1$ als Funktion in den Variablen $p$ und $w$ , so erhält man für die Verstärkungsfaktoren

$k_{11}(p, w) = \left| \frac{p}{p + w} \right| = \underbrace{\left| \frac{1}{1 + (w/p)} \right|}_{\gg 1},$

$k_{12}(p, w) = \left| \frac{w}{p + w} \right| = \underbrace{\left| \frac{1}{1 + (p/w)} \right|}_{\gg 1}.$

Also ist die Variante A im Fall (1.10) nicht stabil. Bei Variante B erhält man dagegen

k 11 (q,λ 2) = k 12 (q,λ 2) = 1.

D. h., der Algorithmus B ist stabil. Für $p : = - 2$ und $q : = 0.01$ ergibt sich bei exakter (bzw. 4-stelliger) Rechnung

$u := p^2 = 4 \quad (4.000),$

$v := u - q = 3.99 \quad (3.990),$

$w := \sqrt{v} = 1.997\, 49\ldots \quad (1.997),$

$\lambda_2 := p - w = -3.997\, 4\ldots \quad (-3.997).$

Die exakte Lösung für $λ 1$ ist $λ 1 = - 0.0025$ . Bei Rechnung mit 4-stelliger Mantisse erhält man im Fall der Variante B $λ 1 = - 0.0025$ , während man für die Variante A $λ 1 = - 0.0030$ erhält mit einem relativen Fehler bezüglich der exakten Lösung von

$\left| \frac{0.0030 - 0.0025}{0.0025} \right| = 0.2$

Im Fall

$p^2 - q > 0, \quad |q| \ll p^2, \quad p > 0$

gilt offenbar $w \approx p$ . Somit ist die Berechnung von $λ 1 : = p + w$ stabil möglich und von $λ 2 : = p - w$ kritisch. Ein stabiler Algorithmus ergibt sich hier durch die Vertauschung der Rollen von $λ 1$ und $λ 2$ oben.

Wir nennen einen Algorithmus zur Lösung eines bestimmten Problems numerisch stabiler als einen zweiten zur Lösung desselben Problems, wenn der Gesamteinfluss aller Rundungsfehler auf die Lösung bei dem ersten Algorithmus kleiner als bei dem zweiten ist.

Kurs:Numerik I/1 Besonderheiten des numerischen Rechnens