\section*{Gleitkomma-Arithmetik} Für $e_{min}, e_{max} \in \mathbb{Z}$, $e_{min} < e_{max}$ ist ein Gleitkommasystem wie folgt definiert: \vspace*{-4mm} \begin{align*} \mathcal{F} &= \mathcal{F}(\beta,t,e_{min},e_{max}) \\ &= \{ \pm m \beta^{e-t} | m \in \N, \beta^{t-1} \leq m \leq \beta^t - 1 \lor m = 0, \\ & \hspace*{16mm}e_{min} \leq e \leq e_{max} \} \end{align*} $x \in \mathcal{F} \setminus \{0\} \Rightarrow \beta^{e_{min}-1} \leq |x| \leq \beta^{e_{max}}(1-\beta^{-1})$. \subsection*{Normalisierte Darstellung} Für $d_1 \neq 0$, $0 < d_1 \leq \beta - 1$: $x=\pm \beta^e ( \frac{d_1}{\beta^1} + \frac{d_2}{\beta^2} + \cdots + \frac{d_t}{\beta^t} ) =: \pm 0.d_1 d_2 \cdots d_t \cdot \beta^e$ \subsection*{Relative Maschinengenauigkeit} $fl(x) \in \mathcal{F}$ ist die $x \in \R$ am nächsten liegende Gleitkommazahl. Für relative Maschinengenauigkeit $\epsilon := \frac{1}{2} \beta^{1-t}$: $\frac{|fl(x)-x|}{|x|} < \epsilon$, $\frac{|fl(x)-x|}{|fl(x)|} \leq \epsilon$ \subsection*{Arithmetische Grundoperationen} Für $x, y \in \mathcal{F}$ sind Operationen $o \in \{x,-,*,\div\}$ bzgl. eines Gleitkommasystems definiert: $\tilde o(x,y) := fl(o(x,y))$ Zu beachten ist hier die Ungültigkeit der Assoziativ- und Distributivgesetze. \subsection*{Kondition mathematischer Probleme} Ein mathematisches Problem $(f,x)$ ist die Auswertung von $f(x)$ an $x \in E$ wobei $f : E \subset X \to R \subset Y$. $X$ und $Y$ sind normierte Räume, $E$ Menge der Eingaben, $R$ Menge der Resultate. Ein Algorithmus für $f$ ist Abbildung $\tilde f : E \subset X \to Y$ s.d. $\tilde f(x)$ in endlich vielen Schritten auswertbar ist und $\tilde f(x) \approx f(x)$ gilt. Die Konditionszahl eines mathematischen Problems $(f, x)$ ist die kleinste Zahl $\kappa_f(x) \geq 0$ mit: \vspace*{-4mm} \[ \frac{\|f(x + \Delta x) - f(x) \|_Y}{\|f(x)\|_Y} \leq \kappa_f(x) \frac{\|\Delta x\|_X}{\|x\|_X} + o(\|\Delta x \|_X) \] Für $\|\Delta x\|_X \rightarrow 0$. \vspace*{-4mm} \[ \kappa_f(x) = \limsup_{\delta \to 0} \left\{ \frac{\|f(x+\Delta x) - f(x)\|_Y \|x\|_X}{\|f(x)\|_Y \|\Delta x\|_X} \right\} \] Für $\Delta x \in X$, $x + \Delta x \in E$, $\|\Delta x\|_X \leq \delta$. Ein Problem $(f, x)$ ist gut konditioniert für \emph{kleine} und schlecht konditioniert für \emph{große} Konditionszahlen $\kappa_f(x)$. Existiert $\limsup$ nicht wird $\kappa_f(x)=\infty$ gesetzt und das Problem als \emph{schlecht gestellt} bezeichnet. \subsubsection*{Kondition stetig differenzierbarer Fkt.} Für $f \in \mathcal{C}^1(E, \R^m)$ in Umgebung $E \subseteq \R^n$ von $x$: \vspace*{-2mm} \[ \kappa_f(x) = \frac{\|f'(x)\|_\infty \cdot \|x\|_X}{\|f(x)\|_Y} \] \subsection*{Stabilität numerischer Algorithmen} Elementaroperationen eines Computers sind mit dem relativen Fehler $\epsilon$ behaftet. Zusätzlich existiert ein Eingabefehler derselben Größenordnung. Es ist also in jedem Fall mit einem relativen Fehler $\kappa_f(x)\epsilon$ zu rechnen. \subsubsection*{Vorwärtsanalyse} \emph{Stabilitätsindikator der Vorwärtsanalyse} eines Algorithmus' $\tilde f$ zur Lösung von $(f,x)$ ist minimales $\sigma = \sigma(x) \geq 0$ für $\{x^\epsilon\}_{\epsilon > 0}$ mit $\|x-x^\epsilon\|_X \leq \epsilon \|x\|_X$: \vspace{-4mm} \[ \frac{\|\tilde f(x^\epsilon) - f(x^\epsilon)\|_Y}{\|f(x^\epsilon)\|_Y} \leq \sigma \kappa_f(x^\epsilon)\epsilon + o(\epsilon) \text{ für } \epsilon \to 0 \] Algorithmus $\tilde f$ ist stabil im Sinne der Vorwärtsanalyse, wenn $\sigma \leq$ \#Elementaroperationen. \subsubsection*{Rückwärtsanalyse} \emph{Stabilitätsindikator der Rückwärtsanalyse} ist minimales $\varrho = \varrho(x) \geq 0$ für $\{x^\epsilon\}_{\epsilon > 0}$ s.d. für bel. $\|x^\epsilon - x\|_X \leq \epsilon \|x\|_X$ Schar $\{\hat x^\epsilon\}_{\epsilon > 0}$ ex. mit $f(\hat x^\epsilon) = \tilde f(x^\epsilon)$: \vspace{-2mm} \[ \frac{\|\hat x^\epsilon - x^\epsilon\|_X}{\|x^\epsilon\|_X} \leq \varrho\epsilon + o(\epsilon) \text{ für } \epsilon \to 0 \] Vorwärtsstabilität folgt aus Rückwärtsstabilität. \section*{Vektor- und Matrixnormen} \subsection*{Induzierte Matrixnorm / Operatornorm} Für Normen $\| \cdot \|_\circ$, $\| \cdot \|_\star$ auf $\K^n$ bzw. $\K^m$ ist eine Matrixnorm $\| \cdot \| : \K^{m \times n} \rightarrow [0,\infty)$ auf dem Vektorraum der $m \times n$-Matrizen definiert: \vspace*{-4mm} \[ \|A\| := \max_{v \in \K^n \setminus \{0\}} \frac{\|Av\|_\star}{\|v\|_\circ} = \max_{\{v \in \K^n | \|v\|_\circ = 1 \}} \|Av\|_\star \] \subsubsection*{Eigenschaften} Für $A \in \K^{m \times n}$ gilt $\forall v \in \K^n : \|Av\|_\star \leq \|A\| \cdot \|v\|_\circ$ Submultiplikativität: $\|AB\| \leq \|A\| \cdot \|B\|$ \subsubsection*{Matrix-$p$-Normen} Induzierte Matrixnorm bei Wahl der $p$-Normen über $\K^n$ bzw. $\K^m$: \vspace*{-4mm} \[ \|A\|_p := \max_{v \in \K^n} \frac{\|Av\|_p}{\|v\|_p} = \max_{\|v\|_p = 1} \|Av\|_p \text{ für } 1 \leq p \leq \infty \] \subsubsection*{Spaltensummennorm} Für $A = (a_1, \cdots, a_n)$ mit $a_j \in \K^m$: \vspace*{-4mm} \[ \|A\|_1 = \max_{1 \leq j \leq n} \|a_j\|_1 = \max_{1 \leq j \leq n} \sum_{i=1}^m |a_{i,j}| \] \subsubsection*{Zeilensummennorm} \vspace*{-4mm} \[ \|A\|_\infty = \max_{1 \leq i \leq m} \sum_{j=1}^n |a_{i,j}| \] \subsubsection*{Spektralnorm} Die Matrix-$2$-Norm wird so genannt, da $\|A\|_2 = \sqrt{\lambda_{max}(A^H A)}$ für $\lambda_{max}(A^H A)$ als Bezeichner des größten Eigenwerts von $A^H A \in \K^{n \times n}$. $\|A\|_2 = \|A^H\|_2$, $\|A^H A\|_2 = \|A\|_2^2$ $\|Q A\|_2 = \|A\|_2$ für unitäre $Q$. \subsection*{Induzierte Normen} Für $A \in \R^{n \times n}$ mit $A > 0$ ist $\skp{z,z}_A := \skp{Az,z}_2$ ein Skalarprodukt auf $\R^n$. Dieses induziert die Energienorm $\|z\|_A = \sqrt{\skp{z,z}_A}$. \subsection*{Kondition einer Matrix} Für $A \in \K^{n \times n} \in GL_n{\R}$, $\|\cdot\|$ induzierte Matrixnorm: \vspace*{-4mm} \begin{align*} \kappa(A) &= \|A\| \cdot \|A^{-1}\| \\ 1 = \|Id\| = \|AA^{-1}\| &\leq \|A\| \cdot \|A^{-1}\| = \kappa(A) \end{align*} \section*{Besondere Matrizen} \subsection*{Diagonaldominante Matrizen} $A \in \R^{n \times n}$ ist diagonaldominant, falls: $\forall i \in \{1,\cdots,n\} : |a_{i,i}| > \sum_{j=1,j\neq i}^n |a_{i,j}|$ Insbesondere sind solche $A$ regulär. Gilt nur $\geq$ so heißt die Matrix schwach diagonaldominant. \subsection*{Positiv definite Matrizen} $A \in \R^{n \times n}$ ist positiv definit d.h. $A > 0$ falls $A=A^T$ und $\forall x \in \R^n \setminus \{0\} : x^TAx > 0$. \subsection*{Hessenberg-Matrizen} Fast obere / untere Dreiecksmatrix wobei 1. untere / obere Nebendiagonale besetzt sein kann. \subsection*{Neumannsche-Reihe} \vspace{-2mm} \[ (Id-M)^{-1} = \sum_{k=0}^\infty M^k \] \subsection*{Bezüglich $A > 0$ konjugierte Vektoren} Vektoren $p, q \in \R^n$ sind konjugiert bzgl. $A > 0$ d.h. $A$-orthogonal, falls $Ap \perp q$, also $\skp{Ap,q}_2=\skp{p,q}_A=0$. \section*{Direkte Verfahren zur LGS Lösung} \subsection*{Cramersche Regel} Sei $A = (a_{i,j})_{ij} \in GL_n(\R)$, $b \in \R^n$, $A[j] = (a_1, \cdots, a_{j-1}, b, a_{j+1}, \cdots, a_n) \in \R^{n \times n}$, $a_k$ k-ter Spaltenvektor von $A$. Dann bildet $x_j = \frac{det(A[j])}{det(A)}$ für $j = 1, \cdots, n$ die eindeutige Lösung $x \in \R^n$ s.d. $Ax=b$. Aufgrund des hohen Aufwands von allg. mehr als $(n+1)!$ arithmetischen Operationen ist die Cramersche Regel nur von theoretischer Bedeutung. \subsection*{Lösung gestaffelter Systeme} Obere Dreicksmatrizen können mittels Rückwärtssubstitution, untere Dreiecksmatrizen mittels Vorwärtssubstitution in $\mathcal{O}(n^2)$ gelöst werden. \subsection*{LR-Zerlegung} $A = LR$ wobei $L$ untere Dreiecksmatrix mit $1$-Diagonale und $R$ obere Dreicksmatrix. \subsubsection*{Berechnung LR-Zerlegung} Die LR-Zerlegung existiert insofern die Diagonaleinträge nicht verschwinden. Insbesondere gilt dies für diagonaldominante Matrizen. \begin{enumerate} \item Spaltenweises nullen der der unteren Einträge mittels \emph{Gauß}, Matrizen $L_1, \cdots, L_{n-1}$ \item $L = L_1^{-1} \cdots L_{n-1}^{-1}$, $R=L_{n-1} \cdots L_1 A$ \end{enumerate} \subsubsection*{Lösung $Ax=b$ mittels LR-Zerlegung} \begin{enumerate} \item $A=LR$ berechnen \item $Lz=b$ Vorwärtssubstitution \item $Rx=z$ Rückwärtssubstitution \end{enumerate} \subsubsection*{Spaltenpivotsuche} Die normale LR-Zerlegung ist nur Vorwärts- und nicht Rückwärtsstabil. Dies kann durch Spaltenpivotsuche verbessert werden. Hier wird in jedem Schritt mittels einer Permutationsmatrix immer mit der größten verbleibenden Zeile eliminiert. \vspace{1mm} Für alle regulären Matrizen existiert eine Spaltenpivotsuchen LR-Zerlegung so, dass $PA=LR$. \subsection*{Cholesky-Zerlegung} Für symmetrische $A > 0$ existiert untere Dreiecksmatrix $L$ mit positiver Diagonale, so dass $A = LL^T$. Sym. $A > 0$ können eindeutig als $A=LDL^T$ geschrieben werden s.d. $L$ untere Dreiecksmatrix mit 1-Diagonale, $D$ positive Diagonalmatrix. $D=D^{1/2}D^{1/2}$ und $G=LD^{1/2}$ ergibt die äquivalente Zerlegung $A=GG^T$. \subsection*{QR-Faktorisierung} Für alle $A \in \R^{m \times n}$ mit $m \geq n$ und $Rang(A)=n$ existiert $A=QR$ mit unitärem $Q \in \R^{m \times m}$ und injektiver oberer Dreiecksmatrix $R$. \subsubsection*{Householder-Reflexionen} \[ H(v) := Id_m - 2 \frac{vv^T}{v^Tv} = Id_m - 2 \frac{vv^T}{\|v\|_2^2} \text{ für } \forall v \in \R^m \setminus \{0\} \] Solche $H(v)$ sind orthogonal und symmetrisch, d.h. $H(v)^T H(v)=Id_m$ und $H(v)^2=Id_m$. Wegen $H(v)v=v-2v=-v$ und $\forall w \in spann\{v\}^\perp : H(v)w=w$ ist $H(v)$ Spiegelung an der Hyperebene $spann\{v\}^\perp$. Solche Reflexionen können durch wiederholte Anwendung Matrizen in obere Dreiecksgestalt überführen: \vspace{1mm} Gesucht ist $v \in \R^m$ für $y \in \R^m$ s.d.: \vspace{-2mm} \[ H(v)y=y - 2 \frac{\skp{v,y}}{\|v\|_2^2}v \overset{!}{=} \alpha e_1 \] Vermeidung von Auslöschung: $\alpha := -sign(y_1)\|y\|_2$ Es ergibt sich mit $v:=y-\alpha e_1$: $H(y-\alpha e_1)y=\alpha e_1$. \vspace{1mm} Seien $Q_k$ die sukzessiven, auf $m \times m$ erweiterten, Householder-Reflexionen. Dann gilt: \vspace{1mm} $R:=Q_p \cdots Q_1 A$, $Q:=Q_1^T \cdots Q_p^T$ s.d. $A=QR$. \subsubsection*{Givens-Rotationen} Mit $c^2 + s^2 = 1, c, s \in \R$ und $l < k$: \[ G(l,k) := \left(\begin{smallmatrix} 1 & & & & & & & & & & \\ & \diagdown & & & & & & & & & \\ & & 1 & & & & & & & & \\ & & & c & & & & s & & & \\ & & & & 1 & & & & & & \\ & & & & & \diagdown & & & & & \\ & & & & & & 1 & & & & \\ & & & -s & & & & c & & & \\ & & & & & & & & 1 & & \\ & & & & & & & & & \diagdown & \\ & & & & & & & & & & 1 \end{smallmatrix}\right) \] Wobei $c$ das Diagonalelement der $l$-ten und $k$-ten Zeile, $s$ $k$-tes Element der $l$-ten Zeile, $-s$ $l$-tes Element der $k$-ten Zeile. Givens-Rotationen sind orthogonal und nicht symmetrisch. $G(l,k)A$ unterscheidet sich von $A$ nur in der $l$-ten und $k$-ten Zeile. \vspace{-4mm} \[ (G(l,k)x)_i = \begin{cases} cx_l + sx_k & i=l \\ -sx_l + cx_k & i=k \\ x_i & \text{sonst} \end{cases} \] $\exists \varphi \in (0,2\pi] : c=\cos{\varphi}, s=\sin{\varphi}$ d.h. $G(l,k)$ ist Rotation um $\varphi$ in Ebene $spann\{e_l,e_k\}$. Ziel: $k$-te Komponente von $x$ nullen für $x_l^2+x_k^2 \neq 0$. \begin{align*} |x_l| > |x_k| : &\tau := \frac{x_k}{x_l}, c := \frac{1}{\sqrt{1+\tau^2}}, s := c\tau \\ |x_l| \leq |x_k| : &\tau := \frac{x_l}{x_k}, s := \frac{1}{\sqrt{1+\tau^2}}, c := s\tau \end{align*} Mit einer solchen Givens-Rotation können einzelne Matrixelemente genullt und $A \in \R^{m \times n}$ so sukzessive in eine obere Dreiecksmatrix transformiert werden. Für Hessenberg-Matrix $A$ ergibt sich $A=QR$ mit: $Q^T := G(n-1,n) \cdots G(1,2)$ und $R:=Q^T A$. \spacing QR mit Householder ist ungefähr doppelt so schnell wie mit Givens. Diese sind daher nur bei strukturierten Matrizen wie Tridiagonal- oder Hessenberg-Matrizen sinnvoll einzusetzen. \section*{Lineare Ausgleichsprobleme} Finde $u^* \in \R^n$: $\|Au^*-b\|_2 = \min_{u\in \R^n} \|Au-b\|_2$ Es sind äquivalent: \begin{enumerate}[label=(\alph*)] \item $u^*$ löst Ausgleichsproblem \item $u^*$ löst Normalengleichung $A^TAu^*=A^Tb$ \item $u^*$ erfüllt $Au^* = P_Ab$ mit Ortogonalprojekt. $P_A : \R^m \rightarrow \R^m$ auf Bild von $A$ \end{enumerate} Das Residuum $Au^*-b$ steht normal zu Bild von $A$. Ein Ausgleichsproblem ist eindeutig lösbar gdw. $Kern(A) = \{0\}$. Die Lösung mit minimaler euklidischer Norm wird Minimum-Norm-Lösung $u^+$ genannt. \subsection*{Singulärwertzerlegung} Sei $A \in \R^{m \times n}$ mit $r=Rang(A) \leq \min\{m,n\}$. Dann existiert $A=USV^T$ mit orthogonalen $U \in \R^{m \times m}$, $V \in \R^{n \times n}$ sowie Diagonalmatrix $S \in \R^{m \times n}$. \section*{Iterative Verfahren zur LGS Lösung} Ein iteratives Verfahren zur Lösung von $Au=b$ liefert zu Startvektor $u^0 \in \R^n$ eine Folge $\{u^k\}_{k\in \N_0} \subset \R^n$ mittels $\Psi_k : (\R^n)^{k+1} \rightarrow \R^n$ durch $u^{k+1} = \Psi_k(u^0, \cdots, u^k)$. Ein Verfahren konvergiert falls $\forall u^0 \in \R^n$, $b \in \R^n : \lim_{k\to \infty} u^k = A^{-1}b$ \subsection*{Abbruchkriterium} Direkte Löser liefern bei exakter Arithmetik nach endlich vielen Schritten $A^{-1}b$. Bei iterativen Lösern gilt im allg. $\forall k \in \N : u^k \neq A^{-1}b$. Ein Abbruchkriterium mit Toleranz $\tau > 0$ ist z.B. $\|u^m-u^{m-1}\| \leq \tau$ oder $\frac{\|Au^m-b\|}{\|b\|} \leq \tau$ (rel. Residuum) \subsection*{Lineare Iterationsverfahren} $u^{k+1} = u^k + N(b-Au^k) = (Id - NA)u^k + Nb$ zu $u^0 \in \R^n$ wobei $N \in GL_n(\R)$ das jeweilige Verfahren charakterisiert. $M := Id - NA$ heißt Iterationsmarix. \subsubsection*{Spektralradius} $\varrho(M) = \max\{|\lambda| : \lambda \in Spec(M)\}$ \subsubsection*{Konvergenz} Ein lineares Iterationsverfahren konvergiert gdw. Spektralradius von $M$ $\varrho(M) < 1$. Ist $\|M\| < 1$ für induzierte Matrixnorm $\|\cdot\|$ dann konvergiert das entsprechende lineare Iterationsverfahren ebenso. Sind $A$ oder $A^T$ diagonaldominant so konvergieren sowohl das Jacobi- als auch das Gauß-Seidel-Verfahren. \subsubsection*{Gesamtschritt- / Jacobi-Verfahren} \[ u_i^{k+1} = u_i^k + \frac{1}{a_{i,i}}\left(b_i - \sum_{j=1}^n a_{i,j} u_j^k \right) \text{ für } i = 1, \cdots, n \] Für $A = D - L - R$: \vspace{-4mm} \begin{align*} u^{k+1} &= u^k + D^{-1}(b-Au^k) \\ &= (Id - D^{-1}A)u^k + D^{-1}b \end{align*} Also: $M^{Jac} = Id - D^{-1}A = D^{-1}(L+R)$, $N^{Jac} = D^{-1}$ \subsubsection*{Einzelschritt- / Gauß-Seidel-Verfahren} \[ u_i^{k+1} = u_i^k + \frac{1}{a_{i,i}}\left(b_i - \sum_{j=1}^{i-1} a_{i,j} u_j^{k+1} - \sum_{j=i}^n a_{i,j} u_j^k \right) \] Für $A = D - L - R$: $M^{GS} = (D - L)^{-1}R$, $N^{GS} = (D - L)^{-1}$. \subsection*{Krylov-Raum-Verfahren} Der Krylov-Raum $m$-ter Ordnung bzgl. $B$ und $v$: \vspace{-4mm} \[ \mathcal{U}_m(B,v) := spann\{v,Bv,B^2v, \cdots, B^{m-1}v\} \subset \R^n \] Das Residuuum der $k$-ten Iterierten $u^k$: \vspace{-2mm} \[ r^k=(I-AN)r^{k-1} = (I-AN)^kr^0 \] Es gilt: $u^k \in V_k$ mit $V_k = u_0 + \mathcal{U}_k(NA,Nr^0)$ sowie $r^0 = b - Au^0$. Minimaleigenschaft der $k$-ten Iterierten bzgl. Norm $\|\cdot\|_\star$ auf $\R^n$: \vspace{-2mm} \[ u^k = argmin\{\|v - A^{-1}b\|_\star : v \in V_k\} \] Ein Krylov-Raum-Verfahren bzgl. einer Norm $\|\cdot\|_\star$ ist nur dann sinnvoll, wenn $u^k$ mit geringem, d.h. im Bereich von zwei Matrix-Vektormultiplikationen liegendem, numerischen Aufwand aus $u^{k-1}$ hervorgeht. \vspace{-2mm} \[ u^{k+1} = u^k + N(b-Au^k) \] $u^k$ wird in jeder Iteration entsprechend der jeweiligen Charakterisierung minimiert. \subsection*{Vorkonditionierer} Anforderungen: $Nv$ sollte schnell zu berechnen sein, weiterhin sollte $N \approx A^{-1}$ möglichst gelten. \subsection*{CG-Verfahren} Das Verfahren der konjugierten Gradienten ist durch die Energienorm charakterisiert und definiert als: \vspace{-2mm} \[ u^k = argmin\{\|v-A^{-1}b\|_A : v \in V_k\} \] Für positiv definite $A, N \in \R^{n \times n}$ sowie $b \in \R^n$ liefert das CG-Verfahren nach spätestens $n$ Schritten $A^{-1}b$. Eigentlich ist es bei exakter Arithmetik also ein direkter Verfahren, wird aber durch früheren Abbruch als iteratives Verfahren verwendet. \subsubsection*{Umsetzung} Minimum $\hat x$ von $f(x)=\frac{1}{2} x^TAx - b^Tx$ ist $A\hat x = b$. $\nabla f(x) = Ax - b \overset{!}{=} 0$ ergibt Minimum da nach oben geöffnete Parabel wg. $A > 0$. Initiale Suchrichtung $p_0$ ist $r_0 = b - Ax_0$. Folgende $p_{k+1}$ ergeben sich aus $p_{k+1} = r_{k+1} - \beta_k p_k$. $p_k$ und $x_k$ sind $A$-konjugiert. $f(x_k + \alpha_k p_k)$ wird via $\alpha_k$ entlang dieser Richtung minimiert mittels: \vspace{-4mm} \begin{align*} \alpha_k &= \frac{\skp{r_k,p_k}_2}{\|p_k\|_A^2} = \frac{\skp{r_k,r_k}_2}{\|p_k\|_A^2} \\ \beta_k &= \frac{\skp{r_{k+1},p_k}_A}{\|p_k\|_A^2} = \frac{\skp{r_{k+1},r_{k+1}}_2}{\|r_k\|_2^2} \end{align*} Es gilt also $x_{k+1} = x_k + \alpha_k p_k$ und $r_{k+1} = r_k - \alpha_k A p_k$. \subsection*{GMRES-Verfahren} Das Verfahren des verallgemeinerten minimalen Residuum liefert die Lösung $Ax=b$ für $A \in GL_n(\R)$ und ist charakterisiert durch: \vspace{-4mm} \begin{align*} u^k :&= argmin\{\|N(b-Av)\|_2 : v \in V_k\} \\ &= argmin\{\|NA(A^{-1}b-v)\|_2 : v \in V_k\} \end{align*} Das GMRES-Verfahren ist also durch die nicht Skalarprodukt-induzierte Norm $\|NA\cdot\|_2$ induziert. \section*{Interpolation} Interpolation von $f$ mit $\varphi$ s.d. $\varphi(t_i) = f(t_i)$ für $i = 0,\cdots, n$. Approximation von $f$ mit $\varphi$ s.d. $\|\varphi - f\|$ möglichst klein in geeigneter Norm. \subsection*{Klassische Polynom-Interpolation} Zu gegebenen Knoten $t_0 < \cdots < t_n$ und Stützwerten $f_i = f(t_i)$ für $i = 0,\cdots,n$ wird Polynom $ \in \Pi_n$ gesucht s.d. $P(t_i)=f_i$ für $i = 0,\cdots,n$. Zu $n+1$ Stützwerten $f_i$ und paarweise verschiedenen Knoten $t_i$ existiert dabei genau ein solches Interpolationspolynom $P=P(f|t_0,\cdots,t_n) \in \Pi_n$. \subsubsection*{Interpolationsfehler} \vspace{-4mm} \[ \|f-P(f|t_0, \cdots, t_n)\|_\infty \leq \sup_{\tau \in [a,b]} \frac{|f^{(n+1)}(\tau)|}{(n+1)!} \|\omega_{n+1}\|_\infty \] $\omega_{n+1} \in \Pi_{n+1}$ ist das \emph{Newton-Polynom} bzgl. $t_0, \cdots, t_n$ mit $\omega_{n+1}(t) := \prod_{i=0}^n (t-t_i)$. \subsection*{Vandermonde-Matrix} \[ \begin{pmatrix} 1 & t_0 & t_0^2 & \cdots & t_0^n \\ 1 & t_1 & t_1^2 & \cdots & t_1^n \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & t_n & t_n^2 & \cdots & t_n^n \end{pmatrix} \begin{pmatrix}a_0 \\ \vdots \\ \vdots \\ a_n\end{pmatrix} = \begin{pmatrix}f_0 \\ \vdots \\ \vdots \\ f_n\end{pmatrix} \] Die Lösung der Vandermonde-Matrix beschreibt $P(f|t_0,\cdots,t_n) \in \Pi_n$, was jedoch zu aufwändig ist. \subsection*{Lagrange-Basis} Basis $\{L_{n,0},\cdots,L_{n,n}\}$ von $\Pi_n$ abhg. $t_0 < \cdots < t_n$ wg. $L_{n,k}(t_i) = \delta_{k,i} = \begin{cases}1 & k=i \\ 0 & sonst\end{cases}$ \[ L_{n,k}(t) := \prod_{j=0,j\neq k}^n \frac{t-t_j}{t_k-t_j} \] Es gilt also: $P(f|t_0,\cdots,t_n)=\sum_{k=0}^n f_k \cdot L_{n,k}$ \spacing Ein Lagrange Polynom zu Stützstelle $t_k$ nimmt an dieser $1$, an allen anderen Stützstellen $0$ an. \subsubsection*{Lemma von Aitken} $P = P(f|t_0,\cdots,t_n)(t) =$ \vspace{-2mm} \[ \frac{(t_0-t)P(f|t_1,\cdots,t_n)(t)-(t_n-t)P(f|t_0,\cdots,t_{n-1})(t)}{t_0-t_n} \] \subsubsection*{Schema von Neville} Sei $t \in \R$ fest, $P_{i,k}(t)=P_{i,k}=P(f|t_{i-k},\cdots,t_i)(t)$ für $i\geq k$. Dann ist insb. $P_{i,0}=f_i$ und $P_{n,n}=P(f|t_0,\cdots,t_n)(t)$ kann rekursiv mit dem \emph{Schema von Neville} berechnet werden: \vspace{-4mm} \begin{align*} P_{i,k} &= \frac{(t_{i-k}-t)P_{i,k-1} - (t_i-t)P_{i-1,k-1}}{t_{i-k}-t_i} \\ &= \frac{t-t_{i-k}}{t_i-t_{i-k}} P_{i,k-1} - \frac{t-t_i}{t_i-t_{i-k}} P_{i-1,k-1} \end{align*} \subsection*{Tschebyscheff-Knoten} Für $i = 0,\cdots,n$: \vspace{-2mm} \[ t_i^{[a,b]} = \frac{b+a}{2} + \frac{b-a}{2} \cos\left(\frac{2i+1}{2n+2} \pi\right) \] Diese Knotenfolge liegt dichter zu den Intervallgrenzen hin und ergibt eine bessere Interpolation als äquidistante Knoten. \subsection*{Satz von Faber} Zu jeder Folge von Knoten $\{t_0^{(n)},\cdots,t_n^{(n)}\}_{n \in \N}$ in $[a,b]$ gibt es ein $f \in C([a,b])$ so, dass $\{P(f|t_0^{(n)},\cdots,t_n^{(n)})\}_{n \in \N}$ für $n \to \infty$ nicht glm. gegen $f$ konvergiert. \section*{Splines} Nachteile der Polynom-Interpolation bei einer größeren Anzahl von Knoten: \begin{enumerate} \item Starke Oszillation des Polynoms \item Konvergenz des Polynoms gegen die interpolierte Funktion ist nicht gewährleistet \end{enumerate} Sei $\Delta = \{t_0,\cdots,t_{l+1}\}$ ein Gitter paarweise verschiedener Knoten $a=t_0 < \cdots < t_{l+1} = b$. $s \in \mathcal{C}^{k-2}(a,b)$ ist Spline der Ordnung $k \in \N$ bzgl. $\Delta$ wenn sie auf jedem Interval $[t_i,t_{i+1}]$ mit einem Polynom $s_i \in \Pi_{k-1}$ übereinstimmt. \subsection*{Spline-Raum} $S_{k,\Delta}$ ist Raum aller Splines der Ordnung $k$ bzgl. $\Delta$. Der Spline-Raum $S_{k,\Delta}$ ist ein reeller Vektorraum mit $\Pi_{k-1}[a,b] \subset S_{k,\Delta}$. Zusätzlich gilt auch $(t-t_i)_+^{k-1} \in S_{k,\Delta}$ \subsubsection*{Abgebrochene Potenzen} Abgebrochene Potenzen vom Grad $k-1$: \[ (t-t_i)_+^{k-1} := \begin{cases}(t-t_i)^{k-1} &: t \geq t_i \\ 0 &: t < t_i\end{cases} \] Für $t_i \in \Delta$, $i \neq l+1$ \subsubsection*{Basis des Spline-Raumes} \[ \mathcal{B} = \{1,t,\cdots,t^{k-1},(t-t_1)_+^{k-1},\cdots,(t-t_l)_+^{k-1}\} \] ist eine Basis von $S_{k,\Delta}$ mit $\dim(S_{k,\Delta}) = k + l$. \subsection*{Spline-Interpolation} Interpolation einer Funktion $f$ bzgl. eines Gitters $\Delta = \{t_0,\cdots,t_{l+1}$ durch Spline der Ordnung $k$. \spacing Im linearen Fall mit $k=2$ stimmt die Anzahl der Knoten $l+2$ mit $\dim(S_{2,\Delta})=l+2$ überein. Es gibt also genau einen Spline der $(t_i,f(t_i))$ interpoliert. \[ I_2 f = \sum_{i=0}^{l+1} f(t_i) B_i \] Für $B_i \in S_{2,\Delta}$ mit $B_i(t_k) = \delta_{i,k}$ \subsection*{Kubische Splines} Kubische Splines der Ordnung 4 eigenen sich für die Darstellung von Kurven, da das menschliche Auge diese als glatt empfindet. \spacing Die Interpolationsbedingungen reichen zur eindeutigen Bestimmung eines interpolierenden Spline aus $S_{4,\Delta}$ nicht aus. Wegen $\dim(S_{4,\Delta}) - (l+2) = l+4-(l+2) = 2$ bleiben zwei Freiheitsgrade unbestimmt. Eine zusätzliche Bedingung ist, dass der interpolierende kubische Spline die minimale Krümmung aller interpolierenden $\mathcal{C}^2$-Funktionen besitzen soll. \subsubsection*{Krümmung einer Funktion} Krümmung von $y : [a,b] \rightarrow \R, y \in \mathcal{C}^2$: \[ \kappa(t) := \frac{y''(t)}{(1+y'(t))^{3/2}} \] $1/\kappa(t)$ ist der Radius des \emph{Krümmungskreises}. Das Krümmungsverhalten von $y$ über ganz $[a,b]$ ist durch ein Integral messbar: \[ \|y''\|_2 := \left(\int_a^b y''(t)^2 dt\right)^{1/2} \] \subsubsection*{Randbedingungen} Für $s \in S_{4,\Delta}$, $\Delta = \{t_0,\cdots,t_{l+1}\}$ sind mögliche Randbedingungen: \begin{enumerate}[label=(\alph*)] \item $s'(a) = f'(a)$ und $s'(b)=f'(b)$: \emph{vollständige Spline-Interpolation} \item $s''(a) = s''(b) = 0$: \emph{natürliche Interpolation} \item $s'(a) = s'(b)$ und $s''(a) = s''(b)$ falls $f$ $b-a$ periodisch: \emph{periodische Interpolation} \end{enumerate} Ist eine dieser Randbedingungen erfüllt, so ist $s$ eindeutig bestimmt. Ferner gilt für alle anderen interpolierenden $y \in \mathcal{C}^2(a,b)$: $\|s''\|_2 < \|y''\|_2$ \subsubsection*{Momente von Splines} Sei $h_{j+1} := t_{j+1} - t_j$ Länge von $[t_j,t_{j+1}]$. $M_j = s''(t_j)$ für $j = 0, \cdots, l+1$ sind die Momente des Splines $s \in S_{4,\Delta}$. Aus den Momenten kann der Spline vollständig rekonstruiert werden. Da $s_j := s|_{[t_j,t_{j+1}]}$ kubisch ist gilt für $s_j'' = s''|_{[t_j,t_{j+1}]}$: \begin{align*} \hspace{-2mm} s_j''(t) &= M_j \frac{t_{j+1}-t}{h_{j+1}} + M_{j+1} \frac{t-t_j}{h_{j+1}} \\ \hspace{-2mm} s_j'(t) &= -M_j \frac{(t_{j+1} - t)^2}{2h_{j+1}} + M_{j+1} \frac{(t-t_j)^2}{2h_{j+1}} + A_j \\ \hspace{-2mm} s_j(t) &= M_j \frac{(t_{j+1} - t)^3}{6h_{j+1}} + M_{j+1} \frac{(t-t_j)^3}{6h_{j+1}} + A_j(t-t_j) + B_j \end{align*} Die Integrationskonstanten $A_j$, $B_j$ lassen sich aus den Interpolationsbedingungen berechnen. \section*{Ein paar Matlab Grundlagen} \begin{lstlisting}[frame=single,language=Matlab] A = [ 1 0 ; 0 1 ] % = eye(2) b = [ 3 4 ]' % = [ 3; 4 ] b(1) = 4 % => b = [ 4 4 ]' A(2,1) = 2 % => A = [ 1 0 ; 2 1 ] c = 1:3 % = [ 1 2 3 ] c = 1:2:6 % = [ 1 3 5 ] A(2,:) % = [ 3 1 ] A(:,1) % = [ 1 ; 3 ] A.^2 % = [ 1 0 ; 9 1 ] \end{lstlisting} \subsection*{Beispiel: LR-Zerlegung} \begin{lstlisting}[frame=single,language=Matlab] function [A, L, R] = lr(A) [w,h] = size(A); if w ~= h error('A is not a square matrix') end for k = 1:w-1 if abs(A(k,k)) < eps error('No LU-decomposition'); end A(k+1:w,k) = A(k+1:w,k) / A(k,k); A(k+1:w,k+1:w) = A(k+1:w,k+1:w) - A(k+1:w,k) * A(k,k+1:w); end L = tril(A,-1) + eye(w); R = triu(A); end \end{lstlisting}