Übungen¶

Aufgabe 1: Gleitkommaarithmetik¶

Welches Ergebnis erhalten Sie für die Rechnung

\(2590 + 4 + 4 + 4\)

in Gleitkommaarithmetik zur Basis \(\beta=10\) mit Mantissenlänge \(t=3\), wenn Sie die Summanden

(a) von links nach rechts und

(b) von rechts nach links

abarbeiten? Welche allgemeine Konsequenz lässt sich daraus ableiten?

Aufgabe 2: Auslöschung¶

Berechnen Sie mit dem Taschenrechner

\(f(x)=\ln(\sqrt{x^2+1}-x)\)

für \(x=10^4\) und \(x=10^6\). Welches Ergebnis ist richtig? Wie können Sie das Ergebnis umformen, um ein richtiges Ergebnis zu erhalten?

Tipp

Vermeiden Sie die Subtraktion von etwa gleich großen Zahlen

Aufgabe 3: Overflow¶

Die Beträge von \(x\) und \(y\) in der Formel

\(r=\sqrt{x^2 + y^2}\)

seien so groß, dass die Auswertung von \(x^2+y^2\) zu einem Overflow führt, obwohl \(r\) im zulässigen Bereich der Maschinenzahlen liegt. Formen Sie die Formel so um, dass ein Overflow vermieden wird.

Aufgabe 4: Zahlentypen¶

a) Stellen Sie fest, mit welchem Zahlentyp Matlab arbeitet, indem Sie realmin, realmax und eps abfragen.

clear;
a = 1;
realmin(a)
realmax(a)
eps(a)

ans =   2.2251e-308
ans =   1.7977e+308
ans =    2.2204e-16

b) Berechnen Sie mit Matlab: \(\sin(0)/\sin(0)\), \(\sin(0)/\cos(\pi/2)\), \(\cos(\pi/2)/\sin(0)\), \(\cos(\pi/2)/\cos(\pi/2)\) und 1+1.0e-18.

c) Versuchen Sie \(eps\) in Matlab selbstständig mit Hilfe einer for- oder while-Schleife zu ermitteln.

Erinnerung: \(eps\) ist die größte Zahl für die gilt: \(x_M=x_M(1+\frac{eps}{2})\), wobei \(x_M\) eine beliebige Maschinenzahl ist.

# calculate eps with a while loop. eps is the biggest number for which x == x(1 + esp/2) is true

Aufgabe 5: Rundungsfehler¶

François Viète hat 1593 als erster eine analytische Gleichung für die Kreiszahl \(\pi\) veröffentlicht. Zu diesem Zeitpunkt gab es zwar bereits einige Ansätze \(\pi\) zu approximieren, doch seine Darstellung als unendliches Produkt

\[\begin{split} \begin{align} a_1 &= \frac{1}{2}\sqrt{2}, \\ a_n &= \frac{1}{2}\sqrt{2+2a_{n-1}}, \text{ für } n \geq 2 \end{align} \end{split}\]

\[ \Rightarrow \frac{2}{\pi} = \lim_{n \to \infty} a_i = a_1 \cdot a_2 \cdot a_3 \cdot ... = \left( \frac{1}{2}\sqrt{2} \right) \cdot \left( \frac{1}{2}\sqrt{2+\sqrt{2}} \right) \cdot \left( \frac{1}{2}\sqrt{2+\sqrt{2 + \sqrt{2}}} \right) \cdot ... \]

war neu. Die Produktformel lässt sich äquivalent als eine Folge schreiben:

\[\begin{split} \begin{align} z_0 &= 2, \notag \\ z_n &= 2^{n+\frac{1}{2}}\sqrt{1-\sqrt{1-4^{-n}z_{n-1}^2}}, \text{ für } n = 1, 2, 3, ... \end{align} \end{split}\]

Die Folge konvergiert für \(n \to \infty\) gegen \(\pi\):

\[ \lim_{n \to \infty} z_n = \pi. \]

Herr Viète konnte damals noch nichts von Computern und den inherent eingebauten Rundungsfehlern wissen. Leider ist gerade seine Folge anfällig für Fehlerfortpflanzung, d.h. kleine Rundungsfehler in jedem Iterationsschritt verstärken sich. Schreiben Sie ein kleines Matlabskript, das diese Fehlerfortpflanzung demonstriert. Berechnen Sie dazu für \(n=1,2,3,...,20\) jeweils die Folgeglieder \(z_n\). Berechnen Sie für jedes \(z_n\) den relativen Fehler

\[ r_n = \frac{|\pi - z_n|}{\pi}. \]

Erstellen Sie einen Plot, in dem der Logarithmus des relativen Fehlers in Abhängigkeit der Iterationszahl \(n\) dargestellt ist. Wie verhält sich der Fehler für kleine bzw. große \(n\)?

% demonstrate numerical error propagation using Viète's formula for approximating pi

Aufgabe 6: Rechengenauigkeit¶

Zunächst schauen wir uns an wie Matlab/Octave Variablen speichert, denen wir einen Wert zu weisen. Dies kann einfach mit dem Befehl \(whos()\) abgefragt werden. Dabei werden alle Variablen, ihre Größen (für Matlab/Octave ist jede Variable eine Matrix), ihr Speicherbedarf sowie ihre Klasse aufgelistet.

clear;
x = 1;
whos

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        x           1x1                          8  double

Total is 1 element using 8 bytes

Die Variable \(x\) ist eine \(1\)x\(1\) Matrix und ist, da wir nichts spezifiziert haben, standardmäßig vom Typ \(double\). Lassen wir uns die Variable ausgeben:

x =  1

wird der Wert jedoch nicht in diesem Datenformat angezeigt. Das Ausgabeformat können wir mit dem Befehl \(format\) ändern:

format longE
x

x =    1.00000000000000E+00

Wie sehen die Ausgaben für andere Formate wie \(shortG\), \(bank\), \(hex\), \(rat\) aus?

Speicherplatz, insbesondere im Arbeitsspeicher ist kostbar und kann bei der Rechnung mit großen Matrizen (n,m >> 1 Mio.) schnell knapp werden und die Performance sinken lassen. Daher kann es Sinn ergeben den Speicherplatz von Variablen zu begrenzen.

Aber Achtung! Wie Sie feststellen werden, kann dies die zuvor untersuchten Fehler verstärken.

Statt double-precision können wir z.B. mit dem Befehl \(single()\) Werte in single-precision abspeichern:

y = single(x);
whos
x
y

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        x           1x1                          8  double
        y           1x1                          4  single

Total is 2 elements using 12 bytes

x =    1.00000000000000E+00
y =    1.00000000000000E+00

In der Ausgabe unterscheiden sich die Werte von \(x\) und \(y\) nicht, der Speicherbedarf von \(y\) ist aber nur halb so groß, wie der von \(x\). Welche Auswirkungen hat das?

Dazu betrachten wir mit dem Befehl \(eps()\) die relative Maschinengenauigkeit der beiden Werte. Der Rundungsfehler kann abgeschätzt werden durch: \(\left|\frac{\operatorname{fl}(x)-x}{x}\right|\leq eps\). Die relative Maschinengenauigkeit kann also interpretiert werden als kleinster Wert, der zu \(x\) addiert und nach der Rundung berücksichtigt wird.

eps(x)
eps(y)

ans =    2.22044604925031E-16
ans =    1.19209289550781E-07

Betrachten wir nun für beide Werte das Phänomen der Auslöschung bei Subtraktion fast gleich großer Zahlen:

z = 1-10^-8 
dx = x-z
dy = y-z

z =    9.99999990000000E-01
dx =    1.00000000502476E-08
dy =    0.00000000000000E+00

Welches Format haben \(dx\) und \(dy\)?

%Show format of dx, dy

Welchen Wert darf \(z\) maximal annehmen, bevor Auslöschung stattfindet? Ab welchem Wert findet Auslöschung auch bei doppelter Genauigkeit statt (\(dx=0\))?

Tipp

Ausprobieren wird Sie viel wertvolle Zeit kosten, haben Sie eine bessere Idee?

Berechnen Sie die relative Maschinengenauigkeit negativer Werte. Was fällt auf? Wie können Sie das erklären?

%Compute eps of negative values

Modellbildung und Simulation

Übungen

Contents