Numerische Mathematik wird genutzt um all jene Probleme der Mathematik zu lösen, die nicht endlich berechenbar sind. Dabei werden die eigentlichen Aufgabenstellungen oft durch geeignete Ersatzprobleme approximiert.

Fehlerbetrachtung

Fast jede Aufgabenstellung der numerischen Mathematik lässt sich als Nullstellenproblem schreiben:

  • gegeben ist eine Funktion $g$ und Eingabedaten $x$
  • suche $y$ mit $g(x,y)=0$

Ist ein Problem gut gestellt, hat es die folgenden Eigenschaften:

  • Existenz: Es gibt mind. eine Lösung
  • Eindeutigkeit: Es gibt zu jedem $x$ genau ein $y$
  • stetige Abhängigkeit: kleine Störungen in $x$ sorgen nur für kleine Änderungen in $y$

Fehlergrößen

  • Aufgabenstellung $x \xrightarrow{f} y$
  • Ersatzproblem $x \xrightarrow{\hat{f}} \hat{y}$
  • Floatingpoint-Implementierung $\tilde{x} \xrightarrow{\tilde{f}} \tilde{y}$
  • $y$ ist exakte Lösung
  • $\tilde{x}$ ist Floatingpoint Darstellung von $x$
  • $\hat{y}$ ist das Ergebnis bei exakter Arithmetix
  • $\tilde{y}$ ist Floatingpointergebnis des Ersatzproblemes mit $\tilde{x}$ als Eingabe
  • Fehler: $e = \tilde{y} - y$
  • absoluter Fehler: $e_a = ||\tilde{y}-y||$
  • relativer Fehler: $e_r = \tfrac{||\tilde{y}-y||}{||y||}$ für $||y|| \neq 0$
  • Schranke für den absoluten Fehler: $\varepsilon _a$
  • Schranke für den relativen Fehler: $\varepsilon _r$
  • Gemischtes Fehlerkriterium: $|| \tilde{y}-y|| \leq \varepsilon _a + ||y|| \varepsilon _r $

Fehlerquellen

Betrachtet wird der Fehler $\tilde{y} -y = \tilde{f}(\tilde{x})-f(x)$.
Umgeformt ergibt sich $\tilde{y} -y = \tilde{f}(\tilde{x})-\hat{f}(\tilde{x})+\hat{f}(\tilde{x})-f(\tilde{x})+f(\tilde{x})-f(x)$

Eingabefehler $f(\tilde{x})-f(x)$

Entsteht bei der Umwandlung der Eingabedaten in Floating-Point-Darstellung $\tilde{x}$. Dieser Fehler tritt immer auf und gilt als unvermeidbar.
Die Größe des Fehlers hängt davon ab, wie empfindlich das Originalproblem $f$ auf Störungen in den Eingabedaten $x$ reagiert.

Verfahrensfehler $\hat{f}(\tilde{x})-f(\tilde{x})$

Hier wird das verwendete Ersatzproblem $\hat{f}$ mit dem Originalproblem $f$ für den selben Input $\tilde{x}$ verglichen.
Die Größe des Fehlers hängt von der Qualität des Ersatzproblems ab.

Akkumulierte Rundungsfehler $\tilde{f}(\tilde{x})-\hat{f}(\tilde{x})$

  • Das Ersatzproblem $\hat{f}$ besteht aus Elementaroperationen
  • Bei der FP-Implementierung des Ersatzproblemes $\tilde{f}$ werden diese Operationen durch ihr FP-Äquivalente ersetzt
  • Die Rundungsfehler der FP-Operationen akkumulieren sich
  • Die Größe des Fehlers hängt von der Robustheit des Ersatzproblemes gegenüber der FP-Implementierung ab.
Floating-Point Darstellung

Festkommadarstellung: Die Anzahl der Nachkommastellen ist fest.

Exponentialdarstellung: Die Zahl wird als Produkt einer Festpunktzahl mit einer Zehnerpotenz dargestellt.

  • Darstellungsform: $v * m * b ^e$
  • Basis: $b$ - natürliche Zahl > 1
  • Vorzeichen: $v$ 1 oder -1
  • Mantisse: $m$ ein $b$-adischer Bruch mit genau einer Vorkommastelle
  • Exponent: $e$ eine ganze Zahl

Es gibt reelle Zahlen, die nicht mit endlicher Mantissenlänge darstellbar sind. (z.B. $\pi$, $e$, oder $\sqrt{2}$)

  • Maschinengenauigkeit: $\varepsilon _M = \tfrac{1}{2} * b ^{-n}$

Das Verhalten von Floating-Point Zahlen auf dem Rechner ist normiert. Es wird die Basis $b = 2$ verwendet und das Verhalten von arithmetischen Grundoperationen ist auch festgelegt.

Single-Format: Vorzeichen + 24-Stellen Mantisse + 8-Stellen Exponent. Etwa 7 gültige Dezimalstellen.

Double-Format: 64-Bit Equivalent zum Single-Format. ungefähr 16 gültige Dezimalstellen.

Eingabefehler, Kondition

Es wird das Verhalten des Eingabefehlers $f(\tilde{x}) - f(x)$ untersucht. Um Zusammenhänge zwischen Fehlern auf der Input-Seite $x$ und der Output-Seite $f(x)$ zu beschreiben.

Die Konditionszahlen beschreiben, wie stark $f$ Störungen in $x$ schlimmstenfalls verstärkt. Bei kleinen Konditionszahlen liegt ein gut konditioniertes Problem vor, bei großen ein schlecht konditioniertes.

Rechenbeispiele und Definitionen folgen

Verfahrensfehler

Das Originalproblem $f$ wird oft erst durch Überführung in ein Ersatzproblem $\tilde{f}$ lösbar. Oft können diverse Ersatzprobleme definiert werden. Der Unterschied zwischen Original- und Ersatzproblem $\hat{f}(\tilde{x})-f(\tilde{x})$ nennt sich Verfahrensfehler. Der Diskretisierungsparameter $h > 0$ beschreibt Qualität und Aufwand des Ersatzproblemes.
Ein vernünftig gewähltes ersatzproblem geht für $h \rightarrow 0$ gegen $0$ (und der Aufwand meist gegen unendlich)

Akkumulierte Rundungsfehler

Wenn unterschiedliche Floating-Point Verfahren genutzt werden kann es zu Ungenauigkeiten (Single- vs. Double-Precision) oder Auslöschungen kommen. Falls Rundungsfehler das Ergebnis stark beeinflussen nennt man die Implementierung instabil, ansonsten stabil.

Rechenbeispiele z.B. zur Auslöschung folgen


  • gut gestellte Probleme besitzen eine eindeutige Lösung, die stetig von den Eingabedaten abhängt
  • schlecht gestellte Probleme werden durch gut gestellte Probleme angenähert
Eingabefehler
  • man versucht, den absoluten oder relativen Gesamtfehler (oder ein gemischtes Kriterium) zu kontrollieren
  • in der Praxis sollten alle drei Fehlerarten in der gleichen Größenordnung gehalten werden, zu große Genauigkeit bei einer Fehlerart allein reduziert den Gesamtfehler nicht nennenswert
    und erhöht eventuell den Aufwand beträchtlich
  • die einzelnen Fehlerarten lassen sich wie folgt charakterisieren
    • Eingabefehler
      • beschreibt die Empfindlichkeit des Originalproblems f gegenüber Störungen in den Eingabedaten x
      • ist bestimmt durch die Kondition von f und damit im Wesentlichen nicht veränderbar
      • gut gestellte aber schlecht konditionierte Probleme können sich genauso verhalten wie schlecht gestellte Probleme
    • Verfahrensfehler
      • beurteilt die Qualität des Ersatzproblems
      • wird in der Regel für jedes Ersatzproblem separat untersucht (z.B. Konvergenztheorie)
      • kann durch ein “besseres” Ersatzproblem kleiner gemacht werden
    • Akkumulierter Rundungsfehler
      • kann durch eine “günstigere” Floating-Point-Implementierung des Ersatzproblems verändert werden (Summationsreihenfolgen, Vermeidung von Auslöschungen etc.)
      • stabile Implementierungen sind robust gegenüber Rundungsfehlereffekten

LGS

Wir betrachten m Gleichungen mit n Unbekannten, die alle nur linear auftreten. Solche Gleichungssysteme können i.A. in der Form $Ax=b$ dargestellt werden. Ist $m=n$ (Gleichungen = Unbekannten) und $A$ regulär, so ist das LGS eindeutig lösbar. A ist regulär, wenn $det(A) \neq 0$, singulär wenn $det(A)=0$.
Die Anzahl unabhängiger Spalten und Zeilen ist immer gleich und wird als $rang(A)$ bezeichnet.

Direkte Löser

I

Gaußelimination
  • Eliminieren durch Skalierbarkeit/Subtraktion von Zeilen unterhalb der Diagonalen
  • erhalte Dreiecksgestalt
  • Lösung durch rückwärts einsetzen

erweitertes System $\xrightarrow {Spaltenelimination} $Dreiecksgestalt $\xrightarrow {Rückwärtseinsetzen}$ Lösung

  • Bei Diagonalelementen = 0 gibt es Probleme, die durch Zeilentausch gelöst werden können
  • Die Gaußelimination ist für jede reguläre Matrix nutzbar.
LU-Zerlegung
Ohne Zeilentausch
  • Gaußelimination zerlegt $A=LU$ unabhäng von rechter Seite $Ax=b$
  • $Ax = b \Leftrightarrow LUx=b \Leftrightarrow Ly=b,Lx=y$
  • $LU$ kann über $A$ gespeichert werden
  • Für mehrere rechte Seiten wird einmal $A=LU$ berechnet und dann wiederholt $LUx=b, LUx'=x'$
Mit Zeilentausch
  • Gaußelimination zerlegt $PA=LU$ unabhäng von rechter Seite $Ax=b$
  • $Ax = b \Leftrightarrow Ly=Pb \Leftrightarrow Ux=y$
  • $P,L,U$ entsehen aus Gaußelimination mit Zeilentausch
  • $L,U$ werden über $A$ gespeichert, $P$ als ganzzahliger Vektor
  • Freiheiten beim Zeilentausch werden später noch ausgenutzt
Cholesky-Zerlegung

Eine Matrix heißt symmetrisch positiv difinit (spd) falls $A=A ^T$ und $< x,Ax > > 0$ gelten. Ist eine Matrix spd so ist sie regulär und alle Diagonalelemente $ d > 0$ .

  • Ist $A$ spd, kann $A$ wie folgt zerlegt werden
    • rationale Cholesky-Zerlegung
      • $A = \tilde{L}D\tilde{L} ^T$
        • $\tilde{L}$ untere Dreiecksmatrix mit 1 auf der Diagonalen
        • $D$ Diagonalmatrix mit positiven Diagonalelementen
    • Cholesky-Zerlegung:
      • $A = LL ^T$
        • $L$ untere Dreiecksmatrix
  • Aufwand Cholesky: ca. $\frac{1}{2}$ Aufwand Gauß, da nur mit symmetrisch Matritzen gearbeitet wird
  • Verschiedene Implementierungen Möglich

II

Fehlerbetrachtung für lineare Gleichungssysteme
  • Betrachte $Ax = b$, $A$ quadratisch und regulär.
  • Problem ist gut gestellt mit $f(A,b)=A ^{-1}b=x$
  • Betrachte zunächst die Kondition : d.h. wie empfindlich reagiert $x$ auf Störungen in $A,b$
  • Das gestörte System sei $\tilde{A}\tilde{x} = \tilde{b}$
  • $\chi(A)=||A||*||A ^{-1}||$ heißt Konditionszahl der Matrix $A$
    • hängt von der verwendeten Norm ab
    • die Äquivalenz der Normen überträgt sich
    • $\chi(A) \geq 1$
    • $\chi(A) = \chi (A ^{-1})$
    • $\chi (AB) \leq \chi (A) * \chi (B)$
  • $\chi(A)$ sollte so klein wie möglich sein
  • Problem: $\chi(A)$ ist in der Regel wegen $||A ^{-1}||$ nicht berechenbar
  • Man sollte vor der Berechnung versuchen, die Kondition zu verbessern
  • $\tilde{A}\tilde{x} = \tilde{b}$ heißt vorkonditioniertes System

  • $\chi (A)$ bestimmt die Konition von $Ax=b$
  • Vorkonditienierung durch Äquilibrieren bzw. approximative Inverse
  • a-posteriori Fehlerkontrolle über Residuum

Pivot-Strategie
  • Die LU-Zerlegung transformiert ein lineares Gleichungssystem in ein anderes
  • $U(x)=y$ einfach lösbar (Rückwärtseinsetzen)
  • wie verändert sich das Problem dadurch ? Wie verändert sich die Kondition (also die Fehlerverstärkung) ? -> Die Elimination erhöht die Fehleranfälligkeit
  • Tausche beim Zeilentausch die Elemente so, dass das betragsgrößte Element auf der Diagonalen steht (Spalten-Pivot-Suche)
  • Effizient, da wenn $A$ spid ist, $g(A) \leq 1$ gilt

  • Benutze Spalten-Pivotsuche um Fehleranfälligkeit zu reduzieren
  • Spalten-Pivotsuche in der Praxis stabil
  • Einfache Implementierung

Orthogonale Transformationen
  • betrachte den k-ten Schritt des LU-Verfahrens: die Spalte ab der Diagonale abwärts wird auf ein Vielfaches des ersten Einheitsvektors transformiert
  • Spaltentransformation auf Vielfaches des Einheitsvektors ist Elementaroperation bei der Matrixzerlegung
  • für orthonormales $Q$ ist $\chi _1 (Q) = 1$ (minimal) und damit $\chi _2 (U) \leq \chi _2 (A)$
  • Drehung und Spiegelung sind orthonormal
  • Im $\mathbb{R} ^2$ können Spalten mit Drehungen eliminiert werden
Givens Verfahren
  • benutze Drehungen zum Eliminieren der Spalten
  • Zu $A$ existiert $Q$ orthonormal mit $A=QU$, $U$ obere Dreiecksmatrix
  • $Q,R$ sind nicht eindeutig, Q i.d.R. voll besetzt

  • Givens benutzt Drehung zum Eliminieren
  • liefert Zerlegung $A=QR$ mit $\chi _2 (R) \leq \chi _2 (A)$
  • Handhabung von Q aufwändiger als bei LU
  • Aufwand ca. 4 * LU
  • Implementierung ähnlich

Householder Verfahren
  • auch Drehungen können zur Elimination genutzt werden
  • Durchführung
    • Betrachte $A=(a _1 ... a _n)$
    • Bestimme im ersten Schritt $Q _1 =I - \tfrac {2}{||v|| ^2 _2} vv ^T , v=a _1 \pm ||a _1|| _2 e _1$ dann gilt $Q _1A= \pm || a _1 || _2e _1$
    • d.h. die erste Spalte ist eliminiert, die Matrix $Q _1$ muss von links auf die restlichen Spalten multipliziert werden
    • wiederholung analog für Schritte $2, ..., n-1$

  • Houshoulder benutzt Spiegelung zur Eliminierung
  • liefert Zerlegung $A=QR$ mit $\chi _2 (R) \leq \chi _2 (A)$
  • Handhabung von $Q$ aufwändiger als bei LU, braucht zusätzlichen Diagonalvektor
  • Aufwand = 2*LU = $\tfrac{1}{2}$ Givens

Iterative Löser

  • Direkte Löser liefern bei exakter Arithmetik nach endlich vielen Schritten die exate Lösung $x$
  • Es gibt auch völlig anders Strukturierte Algorithmen : Iterative Löser
  • Es gibt einstufige und mehrstufige Verfahren. Je nachdem ob einer oder mehrere Vorgänger genutzt werden
  • Wird bei jedem Schritt die selbe Vorschrift genutzt ist das Verfahren stationär, sonst instationär
Lineare Stationäre Verfahren
  • Die Form ist $x _{k+1} = \Phi (x _k)$
  • Es gilt $ x = x _k + A ^{-1} r _k$
    • $A ^{-1}$ ist unbekannt und genau so Aufwendig wie $Ax=b$
    • Idee: Ersetze $A ^{-1}$ durch $M ^{-1}$ mit $M \approx A$ und $M ^{-1}$ einfacher zu berechnen
  • Die Residueniteration mit Vorkonditionierer $M$ ist gegeben durch : $ x _{k+1} = x _k + M ^{-1} r _k$
Jacobi Verfahren
  • Auch Gesamtschrittverfahren genannt
  • Für die Residueniteration gilt : $B = I - M ^{-1} A , d = M ^{-1} b$ bzw. $M = A(I-B)^{-1}$
  • Durchführung
    • $r _k = b - Ax _k$
    • löse $Mp _k = r _k $ d.h. $p _k = M ^{-1} r _k$
    • $x _{k+1} = x _k + p _k$
  • Das Verfahren mit Residueniterations $M=D$ heißt Jacobi Verfahren
  • Als stationäres lineares Iterationsverfahren: $x^{k+1} = Bx^{k+1}+d, B=I-D^{-1}A, d=D^{-1}b$
  • Zerlege $A$ in $A=D-E-F$
  • Sind alle Diagonalelemente ungleich 0, so ist $D$ invertierbar
  • Ist $A$ streng diagonaldominant (Die betragsgrößten Elemente auf der Diagonalen) dann konvergiert das Jacobi Verfahren (und $A$ ist regulär)
  • Ist $A$ spd und $2D-A$ ebenfalls spd so konvergiert das Verfahren sogar monoton ($\rho (B) < 1$)
Gauß-Seidel Verfahren
  • Auch Einzelschrittverfahren
  • Betrachte das Jacobi-Verfahren komponentenweise
  • Der Vorkonditionierer ist hier $M _{GS}=D-E$
  • Sei $A=D-E-F$, $D$ invertierbar, $E+F\geq0$ (komponentenweise) und $\rho (B _j) < 1$ Dann konvergiert das Gauß-Seidel Verfahren schneller als Jacobi.
  • Ist $A$ spd so konvergiert das Verfahren
Nachiteration
  • Je besser $M$ die Matrix $A$ approximiert, desto schneller konvergiert die Residueniteration
  • LU-Verfahren generiert mit Rundungsfehlern $\tilde{L}\tilde{U} \approx A$ und $A\tilde{x}\approx b$
Relaxationsverfahren
  • betrachte das Jacobi-Verfahren $x ^{k+1} = x ^k + D ^{-1}r ^{k}$
  • füge vor dem Korrekturterm einen Parameter $\omega$ ein und Versuche die Konvergenz günstig zu beeinflussen. Das relaxierte Jacobi-Verfahren $J _{\omega} $: $x ^{k+1} = x ^k + \omega D ^{-1}r ^{k}$
  • liegt $\omega$ zwischen 0 und 1 und konvergiert das Jacobi-Verfahren, konvergiert auch das relaxierte
  • $\omega$ wird dabei entweden abgeschätzt oder empirisch bestimmt
Symmetrische Varianten
  • Es gibt Anwendungen (PCG, Multigrid) wo $M$ spd erwünscht ist
  • für das symmetrische Gauß-Seidel-Verfahren (SGS) führt man zunächst einen Gauß-Seidel Schritt und dann einen Rückwärts-Gauß-Seidel Schritt durch
  • Ist $A$ symmetrisch und $a _{ii} > 0$ so ist $M _{SGS}$ spd.

Eigenwerte

Motivation
Theorie
Vektoriteration
Jacobi-Verfahren
QR-Verfahren, Hessenberg-Transformation, Shift-Technik

NLGS

Einfache Verfahren und geometrische Interpolation
Bisektionsverfahren
Regula-Falsi-Verfahren
Sekanten-Verfahren
Taylor-Reihen basierte Verfahren, Newton-Verfahren
Stationäre Verfahren, Bananachscher Fixpunktsatz
Newton Verfahren

Interpolation

Polynominterpolation
Splines
B-Splines
Spline-Kurven und CAD

Approximation

Lineare Ausgleichsprobleme
QR Zerlegung für Lineare Ausgleichsprobleme
CGLS
Pseudoinverse
Singulärwertzerlegung
Regularisierung schlecht konditionierter Probleme
Abgeschnitte Singulärwertzerlegung (TSVD)
Tikhonov Regularisierung
Parameterwahl

Numerische Integration

Newton-Cotes Quadraturformeln
Gauß-Interpolation
Extrapolation
Adaptive Verfahren

Optimierung ohne NB

Abstiegsverfahren
Liniensuche
Steepest-Descent
CG-Verfahren
Newton-Verfahren
Trust-Region Verfahren
Quasi-Newton Verfahren
Nichtlineare Ausgangsprobleme

Diskete Fourier- und Wavelettransformation

Fourierreihe
Diskrete Fouriertransformation
Schnelle Fouriertransformation (FFT)
Diskrete Wavelettransformation