7Ein Beispiel zur Stabilität von Gleitkommaarithmetik
7.1 Berechnung von \(\pi\) nach Archimedes
Eine untere Schranke für \(2\pi\), den Umfang des Einheitskreises, erhält man durch die Summe der Seitenlängen eines dem Einheitskreis eingeschriebenen regelmäßigen \(n\)-Ecks. Die Abbildung links zeigt, wie man beginnend mit einem Viereck der Seitenlänge \(s_4=\sqrt{2}\) die Eckenzahl iterativ verdoppelt.
Abb. 1
Abb.2
Die zweite Abbildung zeigt die Geometrie der Eckenverdoppelung.
Mit \(|\overline{AC}|= s_{n},\quad |\overline{AB}|= |\overline{BC}|= s_{2n},\quad |\overline{MN}| =a, \quad |\overline{NB}| =1-a,\) liefert Pythagoras für die Dreiecke \(MNA\) und \(NBA\) jeweils \[
a^2 + \left(\frac{s_{n}}{2}\right)^2 = 1\qquad\text{bzw.} \qquad
(1-a)^2 + \left(\frac{s_{n}}{2}\right)^2 = s_{2n}^2
\] Elimination von \(a\) liefert die Rekursion \[
s_{2n} = \sqrt{2-\sqrt{4-s_n^2}} \qquad\text{mit Startwert}\qquad
s_4=\sqrt{2}
\] für die Länge \(s_n\)einer Seite des eingeschriebenen regelmäßigen \(n\)-Ecks.
Die Folge \((n\cdot s_n)\) konvergiert monoton von unten gegen den Grenzwert \(2\pi\): \[
n\, s_n \rightarrow 2\pi % \qquad\text{und} \qquad {n c_n}\downarrow 2\pi
\] Der relative Fehler der Approximation von 2π durch ein \(n\)-Eck ist: \[
\epsilon_n = \left| \frac{n\, s_n-2 \pi}{2\pi} \right|
\]
Die Gleichung \[
s_{2n} = \sqrt{2-\sqrt{4-s_n^2}}\qquad \qquad \textrm{(Iteration A)}
\] ist mathematisch äquivalent zu \[
s_{2n} = \frac{s_n}{\sqrt{2+\sqrt{4-s_n^2}}} \qquad \qquad \textrm{(Iteration B)}
\]
(Bitte nachrechnen!)
Allerdings ist Iteration A schlecht konditioniert und numerisch instabil, wie der folgende Code zeigt. Ausgegeben wird die jeweils berechnete Näherung für π.
usingPrintfiterationA(s) =sqrt(2-sqrt(4- s^2))iterationB(s) = s /sqrt(2+sqrt(4- s^2))s_B = s_A =sqrt(2) # Startwert ϵ(x) =abs(x -2π)/2π # rel. Fehler ϵ_A =Float64[] # Vektoren für den Plotϵ_B =Float64[]is =Float64[]@printf""" approx. Wert von π n-Eck Iteration A Iteration B"""for i in3:35push!(is, i) s_A =iterationA(s_A) s_B =iterationB(s_B) doublePi_A =2^i * s_A doublePi_B =2^i * s_Bpush!(ϵ_A, ϵ(doublePi_A))push!(ϵ_B, ϵ(doublePi_B))@printf"%14i %20.15f %20.15f\n"2^i doublePi_A/2 doublePi_B/2end
Während Iteration B sich stabilisiert bei einem innerhalb der Maschinengenauigkeit korrekten Wert für π, wird Iteration A schnell instabil. Ein Plot der relativen Fehler \(\epsilon_i\) bestätigt das:
Weitere Iterationen verbessern das Ergebnis nicht mehr. Sie stabilisieren sich bei einem relativen Fehler von etwa 2.5 Maschinenepsilon:
ϵ_B[end]/eps(Float64)
2.5464790894703255
Die Form Iteration A ist instabil. Bereits bei \(i=16\) beginnt der relative Fehler wieder zu wachsen.
Ursache ist eine typische Auslöschung. Die Seitenlängen \(s_n\) werden sehr schnell klein. Damit ist \(a_n=\sqrt{4-s_n^2}\) nur noch wenig kleiner als 2 und bei der Berechnung von \(s_{2n}=\sqrt{2-a_n}\) tritt ein typischer Auslöschungseffekt auf.
setprecision(80) # precision für BigFloats =sqrt(BigFloat(2))@printf" a = √(4-s^2) als BigFloat und als Float64\n\n"for i =3:44 s =iterationA(s) x =sqrt(4-s^2)if i >20@printf"%i %30.26f %20.16f \n" i x Float64(x)endend
Man sieht die Abnahme der Zahl der signifikanten Ziffern. Man sieht auch, dass eine Verwendung von BigFloat mit einer Mantissenlänge von hier 80 Bit das Einsetzen des Auslöschungseffekts nur etwas hinaussschieben kann.
Gegen instabile Algorithmen helfen in der Regel nur bessere, stabile Algorithmen und nicht genauere Maschinenzahlen!
nach: Christoph Überhuber, „Computer-Numerik“ Bd. 1, Springer 1995, Kap. 2.3↩︎
