JuliaKurs23/chapters/Pi2.qmd

---
engine: julia
---

```{julia}
#| error: false
#| echo: false
#| output: false
using InteractiveUtils
```


# Ein Beispiel zur  Stabilität von Gleitkommaarithmetik



## 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.

:::{.narrow}

| Abb. 1 | Abb.2 |
| :-:    | :-:   |
| ![](../images/pi1.png) | ![](../images/pi2.png)  |
: {tbl-colwidths="[57,43]"}

:::


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|
$$


## Zwei Iterationsvorschriften^[nach: Christoph Überhuber, „Computer-Numerik“ Bd. 1, Springer 1995, Kap. 2.3]
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 π.

```{julia}
using Printf

iterationA(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 in 3:35
    push!(is, i)
    s_A = iterationA(s_A)         
    s_B = iterationB(s_B)         
    doublePi_A = 2^i * s_A
    doublePi_B = 2^i * s_B
    push!(ϵ_A, ϵ(doublePi_A))
    push!(ϵ_B, ϵ(doublePi_B))
    @printf "%14i %20.15f  %20.15f\n" 2^i doublePi_A/2 doublePi_B/2 
end
```

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:


```{julia}
using PlotlyJS

layout = Layout(xaxis_title="Iterationsschritte", yaxis_title="rel. Fehler",
            yaxis_type="log", yaxis_exponentformat="power", 
            xaxis_tick0=2, xaxis_dtick=2)

plot([scatter(x=is, y=ϵ_A, mode="markers+lines", name="Iteration A", yscale=:log10),   
      scatter(x=is, y=ϵ_B, mode="markers+lines", name="Iteration B", yscale=:log10)],
      layout) 
 ```



## Stabilität und Auslöschung

Bei $i=26$ erreicht Iteration B einen relativen Fehler in der Größe des Maschinenepsilons:


```{julia}
ϵ_B[22:28]
```

Weitere Iterationen verbessern das Ergebnis nicht mehr. Sie stabilisieren sich bei einem relativen Fehler von etwa 2.5 Maschinenepsilon:


```{julia}
ϵ_B[end]/eps(Float64)
```


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.


```{julia}
setprecision(80) # precision für BigFloat

s = 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)
    end 
end


```

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!**

:::{.content-hidden unless-format="xxx"}

Offensichtlich tritt bei der Berechnung von $2-a_n$ bereits relativ früh
eine Abnahme der Anzahl der signifikanten Ziffern (Auslöschung) auf, 
bevor schließlich bei der Berechnung von  $a_n=\sqrt{4-s_n^2}$  
selbst schon Auslöschung zu einem unbrauchbaren Ergebnis führt.

:::