JuliaKurs23/nb/Pi2.ipynb

Cell:

[Cell type raw - unsupported, skipped]

Ein Beispiel zur Stabilität von Gleitkommaarithmetik

Berechnung von \pi nach Archimedes

In [ ]:

#| error: false
#| echo: false
#| output: false
using PlotlyJS, Random
using HypertextLiteral
using JSON, UUIDs
using Base64


## see https://github.com/JuliaPlots/PlotlyJS.jl/blob/master/src/PlotlyJS.jl
## https://discourse.julialang.org/t/encode-a-plot-to-base64/27765/3
function IJulia.display_dict(p::PlotlyJS.SyncPlot)
            Dict(
             #   "application/vnd.plotly.v1+json" => JSON.lower(p),
             #   "text/plain" => sprint(show, "text/plain", p),
                "text/html" => let
                    buf = IOBuffer()
                    show(buf, MIME("text/html"), p)
                    #close(buf)
                    #String(resize!(buf.data, buf.size))
                    String(take!(buf))
                end,
                   "image/png" => let
                    buf = IOBuffer()
                    buf64 = Base64EncodePipe(buf)
                    show(buf64, MIME("image/png"), p)
                    close(buf64)
                    #String(resize!(buf.data, buf.size))
                    String(take!(buf))
                end,

            )
        end

   function Base.show(io::IO, mimetype::MIME"text/html", p::PlotlyJS.SyncPlot)
        uuid = string(UUIDs.uuid4())
        show(io,mimetype,@htl("""
    <div style="height: auto" id=\"$(uuid)\"></div>
    <script>
          require(['../js/plotly-latest.min.js'], function(plotly) {        
            plotly.newPlot($(uuid),
            $(HypertextLiteral.JavaScript(json(p.plot.data))),
            $(HypertextLiteral.JavaScript(json(p.plot.layout))),{responsive: true});
            });
        </script>
"""))
end 

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

: {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 π.

In [ ]:

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:

In [ ]:

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:

In [ ]:

ϵ_B[22:28]

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

In [ ]:

ϵ_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.

In [ ]:

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.

:::