234 lines
7.3 KiB
Plaintext
234 lines
7.3 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "0b96655f",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# Ein Beispiel zur Stabilität von Gleitkommaarithmetik\n",
|
||
|
|
"\n",
|
||
|
|
"## Berechnung von $\\pi$ nach Archimedes\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"Eine untere Schranke für $2\\pi$, den Umfang des Einheitskreises, erhält man durch die\n",
|
||
|
|
"Summe der Seitenlängen eines dem Einheitskreis eingeschriebenen regelmäßigen $n$-Ecks.\n",
|
||
|
|
" Die Abbildung links zeigt, wie man beginnend mit einem Viereck der Seitenlänge $s_4=\\sqrt{2}$ die Eckenzahl iterativ verdoppelt.\n",
|
||
|
|
"\n",
|
||
|
|
"| Abb. 1 | Abb.2 |\n",
|
||
|
|
"| :-: | :-: |\n",
|
||
|
|
"|  |  |\n",
|
||
|
|
"\n",
|
||
|
|
"Die zweite Abbildung zeigt die Geometrie der Eckenverdoppelung.\n",
|
||
|
|
"\n",
|
||
|
|
"Mit\n",
|
||
|
|
"$|\\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\n",
|
||
|
|
" $NBA$ jeweils\n",
|
||
|
|
"$$\n",
|
||
|
|
" a^2 + \\left(\\frac{s_{n}}{2}\\right)^2 = 1\\qquad\\text{bzw.} \\qquad\n",
|
||
|
|
" (1-a)^2 + \\left(\\frac{s_{n}}{2}\\right)^2 = s_{2n}^2\n",
|
||
|
|
"$$\n",
|
||
|
|
"Elimination von $a$ liefert die Rekursion\n",
|
||
|
|
"$$\n",
|
||
|
|
" s_{2n} = \\sqrt{2-\\sqrt{4-s_n^2}} \\qquad\\text{mit Startwert}\\qquad\n",
|
||
|
|
" s_4=\\sqrt{2}\n",
|
||
|
|
"$$\n",
|
||
|
|
"für die Länge $s_n$ **einer** Seite des eingeschriebenen regelmäßigen\n",
|
||
|
|
"$n$-Ecks.\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"Die Folge $(n\\cdot s_n)$\n",
|
||
|
|
"konvergiert monoton von unten gegen den\n",
|
||
|
|
"Grenzwert $2\\pi$:\n",
|
||
|
|
"$$\n",
|
||
|
|
" n\\, s_n \\rightarrow 2\\pi % \\qquad\\text{und} \\qquad {n c_n}\\downarrow 2\\pi\n",
|
||
|
|
"$$\n",
|
||
|
|
"Der relative Fehler der Approximation von 2π durch ein $n$-Eck ist:\n",
|
||
|
|
"$$\n",
|
||
|
|
" \\epsilon_n = \\left| \\frac{n\\, s_n-2 \\pi}{2\\pi} \\right|\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"## Zwei Iterationsvorschriften^[nach: Christoph Überhuber, „Computer-Numerik“ Bd. 1, Springer 1995, Kap. 2.3]\n",
|
||
|
|
"Die Gleichung\n",
|
||
|
|
"$$\n",
|
||
|
|
" s_{2n} = \\sqrt{2-\\sqrt{4-s_n^2}}\\qquad \\qquad \\textrm{(Iteration A)}\n",
|
||
|
|
"$$\n",
|
||
|
|
"ist mathematisch äquivalent zu\n",
|
||
|
|
"$$\n",
|
||
|
|
" s_{2n} = \\frac{s_n}{\\sqrt{2+\\sqrt{4-s_n^2}}} \\qquad \\qquad \\textrm{(Iteration B)}\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"(Bitte nachrechnen!) \n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"Allerdings ist Iteration A schlecht konditioniert und numerisch instabil, wie der folgende Code zeigt. Ausgegeben wird die jeweils berechnete Näherung für π.\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "fa3be89a",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"using Printf\n",
|
||
|
|
"\n",
|
||
|
|
"iterationA(s) = sqrt(2 - sqrt(4 - s^2))\n",
|
||
|
|
"iterationB(s) = s / sqrt(2 + sqrt(4 - s^2))\n",
|
||
|
|
"\n",
|
||
|
|
"s_B = s_A = sqrt(2) # Startwert \n",
|
||
|
|
"\n",
|
||
|
|
"ϵ(x) = abs(x - 2π)/2π # rel. Fehler \n",
|
||
|
|
"\n",
|
||
|
|
"ϵ_A = Float64[] # Vektoren für den Plot\n",
|
||
|
|
"ϵ_B = Float64[]\n",
|
||
|
|
"is = Float64[]\n",
|
||
|
|
"\n",
|
||
|
|
"@printf \"\"\" approx. Wert von π\n",
|
||
|
|
" n-Eck Iteration A Iteration B\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"for i in 3:35\n",
|
||
|
|
" push!(is, i)\n",
|
||
|
|
" s_A = iterationA(s_A) \n",
|
||
|
|
" s_B = iterationB(s_B) \n",
|
||
|
|
" doublePi_A = 2^i * s_A\n",
|
||
|
|
" doublePi_B = 2^i * s_B\n",
|
||
|
|
" push!(ϵ_A, ϵ(doublePi_A))\n",
|
||
|
|
" push!(ϵ_B, ϵ(doublePi_B))\n",
|
||
|
|
" @printf \"%14i %20.15f %20.15f\\n\" 2^i doublePi_A/2 doublePi_B/2 \n",
|
||
|
|
"end"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "7ce50a85",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"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:\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "5004b70e",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"using CairoMakie\n",
|
||
|
|
"\n",
|
||
|
|
"F = Figure(resolution=(800,500))\n",
|
||
|
|
"ax = Axis(F[1, 1], xlabel = \"Iterationsschritte\", ylabel = \"rel. Fehler\",\n",
|
||
|
|
" yscale = log10, xticks=2:2:35, yticks=LogTicks(-20:0))\n",
|
||
|
|
"\n",
|
||
|
|
"scatterlines!(ax, is, ϵ_A, label = \"Iteration A\" )\n",
|
||
|
|
"scatterlines!(ax, is, ϵ_B, marker=:xcross, label = \"Iteration B\" )\n",
|
||
|
|
"\n",
|
||
|
|
"axislegend(position=:rc)\n",
|
||
|
|
"F"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "7a93c5e3",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Stabilität und Auslöschung\n",
|
||
|
|
"\n",
|
||
|
|
"Bei $i=26$ erreicht Iteration B einen relativen Fehler in der Größe des Maschinenepsilons:\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "527eee0e",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"ϵ_B[22:28]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "6f2f4475",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Weitere Iterationen verbessern das Ergebnis nicht mehr. Sie stabilisieren sich bei einem relativen Fehler von etwa 2.5 Maschinenepsilon:\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "3b2a8b77",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"ϵ_B[end]/eps(Float64)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "6ff2d4b8",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Die Form Iteration A ist instabil. Bereits bei $i=16$ beginnt der relative Fehler wieder zu wachsen.\n",
|
||
|
|
"\n",
|
||
|
|
"Ursache ist eine typische Auslöschung. Die Seitenlängen $s_n$ werden sehr schnell klein. Damit ist\n",
|
||
|
|
"$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.\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "ca6c5ca2",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"setprecision(80) # precision für BigFloat\n",
|
||
|
|
"\n",
|
||
|
|
"s = sqrt(BigFloat(2))\n",
|
||
|
|
"\n",
|
||
|
|
"@printf \" a = √(4-s^2) als BigFloat und als Float64\\n\\n\"\n",
|
||
|
|
"\n",
|
||
|
|
"for i = 3:44\n",
|
||
|
|
" s = iterationA(s)\n",
|
||
|
|
" x = sqrt(4-s^2)\n",
|
||
|
|
" if i > 20\n",
|
||
|
|
" @printf \"%i %30.26f %20.16f \\n\" i x Float64(x)\n",
|
||
|
|
" end \n",
|
||
|
|
"end\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "5e4994b3",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"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. \n",
|
||
|
|
"\n",
|
||
|
|
"**Gegen instabile Algorithmen helfen in der Regel nur bessere, stabile Algorithmen und nicht genauere Maschinenzahlen!**\n",
|
||
|
|
"\n",
|
||
|
|
":::{.content-hidden unless-format=\"xxx\"}\n",
|
||
|
|
"\n",
|
||
|
|
"Offensichtlich tritt bei der Berechnung von $2-a_n$ bereits relativ früh\n",
|
||
|
|
"eine Abnahme der Anzahl der signifikanten Ziffern (Auslöschung) auf, \n",
|
||
|
|
"bevor schließlich bei der Berechnung von $a_n=\\sqrt{4-s_n^2}$ \n",
|
||
|
|
"selbst schon Auslöschung zu einem unbrauchbaren Ergebnis führt.\n",
|
||
|
|
"\n",
|
||
|
|
":::"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Julia 1.8.5",
|
||
|
|
"language": "julia",
|
||
|
|
"name": "julia-1.8"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|