first commit cont.
This commit is contained in:
225
chapters/Pi2.qmd
Normal file
225
chapters/Pi2.qmd
Normal file
@@ -0,0 +1,225 @@
|
||||
---
|
||||
format:
|
||||
html:
|
||||
include-in-header:
|
||||
text: |
|
||||
<script type="application/javascript">
|
||||
window.PlotlyConfig = {MathJaxConfig: 'local'}
|
||||
</script>
|
||||
---
|
||||
|
||||
# Ein Beispiel zur Stabilität von Gleitkommaarithmetik
|
||||
|
||||
## Berechnung von $\pi$ nach Archimedes
|
||||
|
||||
```{julia}
|
||||
#| 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 π.
|
||||
|
||||
```{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.
|
||||
|
||||
:::
|
||||
Reference in New Issue
Block a user