cm.tex

%Header
\documentclass{scrartcl}
\usepackage[ngerman]{babel}
\usepackage[utf8]{inputenc}	%Umlaute
\usepackage{paralist}
\usepackage{amsmath}        % align.
\usepackage[justification=centering]{caption} % Für \caption
\usepackage{graphicx}
\usepackage{chngcntr}
\usepackage{units}
\usepackage{amssymb}
\usepackage[locale = DE]{siunitx}
\usepackage{rotating}
\counterwithin{equation}{section}
\usepackage{listings}             % Include the listings-package
\lstset{language=Pascal}   

\begin{document}
\title{Computergestüzte Methoden der exakten Naturwissenschaften}           
\maketitle
\tableofcontents
\newpage
%Header

\section{Fehler}
Ein Ziel der Naturwissenschaften ist die Beschreibung der Natur mit Hilfe von mathematischen Gleichungen und deren Lösungen, daraus ergibt sich allerdings ein Problem.
\begin{itemize}
\item[\textbf{Problem:}] Die Gleichungen der naturwissenschaftlichen  Beschreibungen können nicht immer mit Bleistift und Papier zu gelöst werden.
\item[\textbf{Lösung 1:}] Vereinfachung der Gleichungen $\hat{=}$ Näherung/Approximation
\item[\textbf{Lösung 2:}] Numerische Lösung der Gleichungen.
\end{itemize}
Diese Vorlesung möchte sich mit der zweiten Lösungsmethode befassen, hierbei ist es allerdings wichtig die Genauigkeit der numerisch ermittelten Ergebnisse (die Fehler) mit zu berücksichtigen.

  Allgemein gibt es für es verschiedene Quellen für Fehler:
\begin{itemize}
 \item[\textbf{Eingabefehler:}] Diese entstehen durch Ungenauigkeiten innerhalb der Eingabedaten.
  \item[\textbf{Näherungsfehler:}] Solche entstehen aus der Verwendung vereinfachter mathematischer Ausdrücke anstelle der exakten.
  \item[\textbf{Modellfehler:}] Diese entstehen aus der Nutzung vereinfachter physikalischer Modelle.
  \item[\textbf{Rundungsfehler:}] Solche entstehen aus der numerischen Darstellung von Zahlen und der damit verbundenen endlichen Genauigkeit.
\end{itemize}

\subsection{Beispiele für Näherungsfehler}
Viele mathematische Gleichungen der Physik sind in ihren exakten Formulierungen nicht oder nur sehr aufwendig lösbar. Ein Ausweg stellen Approximationen dar aus welchen allerdings zusätzliche Näherungsfehler resultieren.
Beispiele hierfür sind über unendliche Reihen definierte Funktionen aber auch Differentialgleichungen im Kontinuum.
\begin{itemize}
\item[\textbf{Exponentialfunktion:}]Die Exponentialfunktion ist definiert durch:
\begin{align*}
e^x=\sum_{n=0}^{\infty}\frac{x^n}{n!}.
\end{align*}
Eine solche Funktion kann durch eine endliche Reihe genähert werden:
\begin{align*}
e^x=\sum_{n=0}^{N}\frac{x^n}{n!}
\end{align*}
\item[\textbf{Differentialgleichung im Kontinuum:}]Eine Differentialgleichung im Kontinuum kann durch die Lösung der zugehörigen diskretisierten Gleichung genähert werden.
Sei die Differentialgleichung gegeben durch:
\begin{align*}
\frac{d}{dx}f(x)=a \ f(x),
\end{align*}
so ergibt sich die dikretisierte Gleichung aus der Diskretisierung auf  
bestimmte Gitterpunkte $x_i$ mit dem Abstand $\Delta x=x_{i+1}-x_i$:
\begin{align*}
\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}=a \ \frac{f(x_{i+1})+f(x_i)}{2}.
\end{align*}
Zur Verbesserung der Diskretisierung kann dann $\Delta x$ immer weiter gegen $0$ gesetzt werden.
Ein \textit{Nachteil} ist hierbei die Erhöhung der Rechenoperationen und der damit verbundenen Rechenzeit. Außerdem vergrößern sich hiermit die Rundungsfehler.
\end{itemize}
Das Ziel der Numerik besteht nun im optimalen Kompromiss zwischen Fehler und Rechenzeit.

\subsection{Beispiel für Modellfehler}
Als Beispiel für einen aus einem Modell resultierenden Fehler wird die Planetenbewegung betrachtet. Nach dem ersten newtonschen Gesetz gilt:
\begin{align*}
\textbf{F}=m \textbf{a}=m \boldsymbol{\ddot{r}}=-\frac{M \ \textbf{r}}{|\textbf{r}|^3}
\end{align*}
Hierbei gehen allerdings eine reihe von Näherungen ein:
\begin{itemize}
\item Die Sonnenmasse $M$ wird relativ zur Planetenmasse als sehr groß angenommen
\item Eine geschwindigkeitsabhängige Reibungskraft $\boldsymbol{F_R}=\gamma \ \boldsymbol{\dot{r}}$ wird vernachlässigt, diese ist allerdings für kleinere Objekte wichtig.
\item Auch relativistische Effekte werden vernachlässigt, solche erklären allerdings Phänomene wie die Periheldrehung des Merkurs.
\item Eigentlich handelt es sich um ein Mehrkörperproblem der Form:
\begin{align*}
m_1 \boldsymbol{\ddot{r_i}}=\sum \boldsymbol F_i=- \sum_{j\neq i} G m_i m_j \frac{\boldsymbol{r_i}-\boldsymbol{r_j}}{|\boldsymbol{r_i}-\boldsymbol{r_j}|^3}.
\end{align*}
Eine Gleichung für mehr als 3 Objekte kann so also leicht geschrieben werden.
Allerdings ist das Problem bereits ab einer Beteiligung von 3 Objekten nur noch unter Annahme bestimmter Bedingungen und ab 4 Objekten überhaupt nicht mehr exakt lösbar.
\end{itemize}

\subsection{Rundungsfehler}
Beim durchführen von Rechenoperationen mit reellen Zahlen am Computer muss gerundet werden, die daraus entstehenden Fehler heißen Rundungsfehler. 
\subsubsection{Gleitpunktarithmetik}
Reelle Zahlen werden am Computer in das Gleitpunktformat umgewandelt.
Der Vorteil gegenüber dem Festpunktformat liegt im geringeren Speicherbedarf.
Hierzu werden die eingegebenen Zahlen in der Form:
\begin{align*}
x=\pm \sum_{i=1}^n z_i \ B^{E-i} := \pm (\underbrace{0,z_1 z_2 ... z_n}_\text{Mantisse})_B \ B^E 
\end{align*}
dargestellt.
Dabei gilt des weiteren für den Exponenten $E \in \mathbb{Z}$: $m \le E \le M$.
Außerdem gilt $z_i \in \{0,1...B-1\}$. \\
\textbf{Beispiel}:\\
\begin{align*}
1234,567=(0,1234567)_{10} \cdot 10^4
\end{align*}

Die Werte für $n$, $B$, $m$ und $M$ sind hierbei maschinenabhängig, werden also durch den Rechner und den Compiler bestimmt. \\
Übliche Basen sind:
\begin{itemize}
\item[$B=2$:] Dualzahlen
\item[$B=8$:] Oktalzahlen
\item[$B=10$:] Dezimalzahlen
\item[$B=16$:] Hexadezimalzahlen
\end{itemize}
\textbf{Standartformate für B=2:}\\
\textbf{Single:}
Dieses Format besteht aus 32 Bits bzw. 4 Bytes.
Diese ergeben sich aus:
\begin{itemize}
\item[Vorzeichen:] 1 Bit
\item[Exponent:] 8 Bits
\item[Mantisse:] 23 Bits
\item[Genauigkeit:] 6 Ziffern unterscheidbar
\end{itemize}
\textbf{Double:}
Dieses Format besteht hingegen aus 64 Bits:
\begin{itemize}
\item[Vorzeichen:] 1 Bit
\item[Exponent:] 11 Bits
\item[Mantisse:] 52 Bits
\item[Genauigkeit:] 15 Ziffern unterscheidbar
\end{itemize}
\textbf{Beispiel}
Binäre Darstellung von $(5,0625)_{10}$:
\begin{align*}
(5,0625)_{10}=(0,50625)_{10} \cdot 10^1=2^2+2^0+2^{-3}+2^{-4}\\
\Rightarrow (5,0625)_{10}=(101,0001)_2=(0,1010001)_2 \cdot 2^{(11)_2}
\end{align*}
Manche Zahlen wie $(0,3)_{10}$ lassen sich allerdings nur schwer als duale Zahlen darstellen.

Die \textbf{größte darstellbare Zahl} ergibt sich zu:
\begin{align*}
x_{max}=(0,\underbrace{[B-1] [B-1] ... [B-1]}_\text{n Ziffern})_B \ B^M=B^M[B-1]\frac{B^{-n}(B^n-1)}{B-1}=B^M(1-B^{-n}).
\end{align*}
Dagegen ergibt sich die \textbf{kleinstmögliche Zahl} zu:
\begin{align*}
x_{min}=B^{m-1}
\end{align*}
Folglich ist die Menge der darstellbaren Maschinenzahlen endlich.
Ergibt sich während der Rechnung eine Zahl $x>x_{max}$ folgt ein overflow und die Zahl wird auf $\infty$ gesetzt.
In gleicher weise ergibt sich für $x<x_{min}$ der underflow und die Zahl wird auf 0 gesetzt.
\textbf{Beispiele:}
\begin{align}
x_{max}+x_{max}=\infty \\
x_{min} \ B^{-1}=0
\end{align}
Jede reelle Zahl die keine Maschinenzahl ist muss in eine solche umgewandelt werden. Idealerweise wählt man Maschinenzahl dabei möglichst nahe der reellen Zahl $\hat{=}$ Rundung.
\subsubsection{Rundung}
Beim Runden wird für eine Zahl $x$ eine Näherung $rd(x)$ unter den Maschinenzahlen geliefert, so dass der absolute Fehler $|x-rd(x)|$ minimal ist. Der dabei unvermeidbare Fehler heißt Rundungsfehler.
Eine n-stellige Dezimalzahl im Gleitpunktformat $\tilde{x}=\pm(0,z_1 ... z_n)_{10}=rd(x)$ hat einen maximalen Fehler von:
\begin{align*}
|x-rd(x)| \leq 0,\underbrace{00...00}_\text{n Ziffern}5 \cdot 10^E=0,5 \cdot 10^{E-n}.
\end{align*}
Für eine allgemeine Basis B ergibt sich:
\begin{align*}
|x-rd(x)|\leq \frac{B}{2} \ \frac{1}{B}  \ B^{E-n}=\frac{1}{2}B^{E-n}.
\end{align*}
Rundungsfehler werden durch die gesamte Rechnung getragen.

Bei einer \textbf{n-stellige Gleitpunktarithmetik}
wird jede einzelne Rechenoperation auf $n+1$ Stellen genau berechnet und dann auf $n$ Stellen gerundet. Es wird also nicht nur das Endergebnis gerundet.\\
\textbf{Beispiel:} $2590+4+4$ in \textit{3-stelliger dezimaler Gleitpunktarithmetik}\\
Von links nach rechts:
\begin{align*}
2590+4=2594\underbrace{\Rightarrow}_\text{Rundung}2590\\
2590+4=2594\underbrace{\Rightarrow}_\text{Rundung}2590\\
\end{align*}
Von rechts nach links:
\begin{align*}
4+4=8\underbrace{\Rightarrow}_\text{Rundung}8\\
2590+8->2598\underbrace{\Rightarrow}_\text{Rundung}2600
\end{align*}
Das exakte Ergebnis wäre $2598$. 
Die Reihenfolge der Ausführungen der Rechenoperationen verändert also das Ergebnis.
Daraus folgt die \textbf{Regel}, dass beim \textbf{Addieren} die Summanden in der Reihenfolge ihrer aufsteigender Beträge addiert werden.
So erhält man bei gleicher Rechenzeit bessere Ergebnisse.\\
\\
\textit{Einschub:Maß für die Rechenzeit eines Computers:\\
$\text{flops} \hat{=} \text{floating point operations per second}$, dabei sind Multiplikation und Division typische Operationen.
Eine Rangliste schnellsten Computer wird auf www.top500.org geführt.}\\
\\
Der \textbf{relative Fehler} ist meist relevanter als der absolute Fehler.
Die Näherung $\tilde{x}$ zu dem exaktem Wert $x$ ergibt einen relativer Fehler: $\epsilon = |\frac{\tilde{x}-x}{x}| \approx |\frac{\tilde{x}-x}{\tilde{x}}|$.\\
Daraus ergibt sich der maximaler Rundungsfehler zu: 
\begin{align*}
\epsilon_{max}=\frac{\frac{1}{2} \ B^{E-n}}{B^{E-1}}=\frac{1}{2} \ B^{1-n}
\end{align*}
%
Für duale Rechnungen im Computer gilt also $B=2 \epsilon_{max} \cdot 2^{-n}$.\\
$\epsilon_{max}$ wird auch \textbf{Maschinengenauigkeit} genannt und gibt die kleinste positive Zahl an für die gilt $1 \cdot \epsilon_{max} \neq 1$.\\
$\epsilon_{max}$ kann aus Rechenoperationen rekonstruiert werden.\\ 
\\
\textbf{Rundungsfehler bei Rechenoperationen}\\
\textbf{Beispiele}: (mit 4er Mantissen und 1er Exponentenziffer, Dez.)\\
\\
\textit{Addition und Subtraktion von Zahlen mit stark unterschiedlichen Exponenten}:
\begin{itemize}
\item[Rundungsfehler kann verloren gehen:]$1234+0,5=\SI{0,1234e4}+\SI{0,5e0}=1235$.\\
Fehler: $0,5$, rel. Fehler: $0,00040$, $\epsilon_{max}=\SI{0,5e-3}+\SI{0,5e0}=1235$, der Rundungsfehler ist also kleiner als der Maximale.
\end{itemize}
\textit{Multiplikation und Division:}
\begin{itemize}
\item[(underflow/overflow möglich):]
$0,2 \cdot 10^{2} \cdot \SI{0,3e-6}=\SI{0,6e-12}=0$\\
$0,2 \cdot 10^{-2} \cdot \SI{0,3e-6}=\SI{0,6e12}=\infty$ (Hier wäre der rel. Fehler $\infty$)
\end{itemize}
\textit{Fehler beim Assoziativgesetz:}
\begin{itemize}
\item[a)] $\SI{0,1111e-3}+(-0,1234+0,1234)=\SI{9,111e-3}+0,0009=\SI{0,10111e-2}=\SI{0,1011e-2})$
\item[b)]$(-0,1234+0,1234= \SI{9,111e-3}+0,0009)+0,1243=-0,1233+0,1243=0,0010=0,100 \cdot 10^{-2}
$\end{itemize}
Der exakte Wert wäre aber: $0,10111 \cdot 10^{-2}$
Daraus folgt:
\begin{itemize}
\item[a)]Fehler: $0,00001 \cdot 10^{-2}$, rel. Fehler: $0,01\%$
\item[b)]Fehler: $0,00111 \cdot 10^{-2}$, rel. Fehler: $1 \%$
\end{itemize}
Im Fall b) ist also $\epsilon>\epsilon_{max}$.


\subsubsection{Fehlerfortpflanzung bei Rechenoperationen}
Fehler werden beim Rechnen weitergetragen, selten werden dabei die Fehler kleiner (meistens werden sie größer!).
Durch das Umstellen von Formeln können Fehler minimiert werden, trotzdem müssen Fehler abgeschätzt werden. \\
\\
\textbf{Additionsfehler:}\\
Gegeben: Fehlerbehaftete Größen $\tilde{x}$ und $\tilde{y}$ zu den Werten $x$ und $y$.\\
Fehler der Summe: $\tilde{x}+\tilde{y}-(x+y)=(\tilde{x}-x)+(\tilde{y}-y)$\\
Im ungünstigsten Fall addieren sich die Fehler: bei Additionen und Subtraktionen addieren sich die Absolutbeträge der Fehler der einzelnen Terme.\\
\\
\textbf{Multiplikation:}
Fehler: $\tilde{x} \tilde{y}- x y=\tilde{x}(\tilde{y}-y)+\tilde{y}(\tilde{x}-x)-(\tilde{x}-x)(\tilde{y}-y)$,
also hat das Produkt von $\tilde{y}$ mit einer Maschinenzahl ohne Fehler den $\tilde{x}$-fachen Fehler; Produkt der Fehler typischerweise vernachlässigbar.\\
Der absolute Fehler eines Produkts ist gegeben durch das Produkt des Faktors mit dem Fehler des anderen Faktors. (=2 Terme, oft ist einer der Fehler dominant.)\\
\\
\textbf{Relativer Fehler Multiplikation:}\\
\begin{align*}
\frac{\tilde{x} \tilde{y}- x y}{\tilde{x \tilde{y}}}=\frac{\tilde{y}-y}{\tilde{y}}+\frac{\tilde{x}-x}{\tilde{x}}-\frac{(\tilde{x}-x)(\tilde{y}-y)}{\tilde{x} \tilde{y}},
\end{align*} beim Multiplizieren addieren sich die relativen Fehler, Division analog.

\subsubsection{Fehlerfortpflanzung bei Funktionen}
Die Funktion wird an der Stelle $\tilde{x}$ anstatt $x$ ausgewertet, daraus folgt ein fehlerbehafteter Funktionswert.
Je nach Funktion resultiert ein kleiner oder großer Fehler.
Bei weiteren Funktionsauswertungen wird der Fehler typischerweise größer.\\
Aus dem Mittelwertsatz folgt:
\begin{align*}
\int_x^{\tilde{x}}g(x')dx'=g(x_0)(\tilde{x}-x)\\
\frac{\int_x^{\tilde{x}}g(x')dx'}{(\tilde{x}-x)}=g(x_0),
\end{align*}
an einer unbekannten Stelle $x_0$ im Intervall $(x,\tilde{x})$.\\
Wähle $g(x)=f'(x)$:
\begin{align*}
|f(\tilde{x})-f(x))|=|\tilde{x}-x||f'(x_0)|.
\end{align*}
Der absolute Fehler vergrößert sich also für $|f'(x_0)|>1$ und wird für $|f'(x_0)|<1$ kleiner.
Die Ableitung kann also als Verstärkungsfaktor des Fehlers interpretiert werden.\\
\\

\textbf{Abschätzung des absoluten Fehlers:}\\
\begin{align*}
|f(x)-f(\tilde{x})|\leq M \ |x-\tilde{x}|,
\end{align*}
mit $M=\underset{x\leq x_0 \leq \tilde{x}} {max} (|f'(x_0)|)$.
Schätzung des Fehlers: $|f(x)-f(\tilde{x})| \approx f'(\tilde{x})|x-\tilde{x}|$.\\
\\
\textbf{Beispiel 1:}\\
Fortpflanzung des absoluten Fehlers von $f(x)=sin(x)$:
\begin{align*}
f'(x)=cos(x) \rightarrow M=1,
\end{align*}
das heißt für die meisten Argumente verringert sich der absolute Fehler.\\
\textbf{Beispiel 2:}
\begin{align*}
f(x)=\sqrt{x};f'(x)=\frac{0,5}{\sqrt{x}},
\end{align*}
divergiert also für $x \rightarrow 0$\\
\\

\textbf{Der relative Fehler bei Funktionsauswertung:}
\begin{align*}
\frac{|f(x)-f(\tilde{x}|)}{|f(x)|}\leq \frac{M |x|}{|f(x)|} \frac{|x-\tilde{x}|}{|x|}\\
\approx \underbrace{\frac{|f'(\tilde{x})||\tilde{x}|}{|f(\tilde{x})|}}_{\textbf{Konditionszahl}} \cdot  \frac{|x-\tilde{x}|}{|\tilde{x}|}.
\end{align*}
Die Konditionszahl ist also ein Verstärkungsfaktor für relative Fehler;
qualitativ: Probleme wenn Konditionszahl$>>$1 ;schlecht konditioniertes Problem.
\section{Nullstellenproblem}
Gegeben: Funktion $\mathbb{R} \rightarrow \mathbb{R}$\\
Gesucht: Nullstellen also $x_0$ aus R mit $f(x_0)=0$
Grundsätzlich:
\begin{itemize}
\item Gibt es überhaupt Nullstellen? Wenn ja, in welchem Bereich?
\item Gibt es mehrere Lösungen?
\end{itemize}

\textbf{Zwischenwertsatz}\\
$f:[a,b]\rightarrow R$, stetig für $C \in R$ mit $f(a) \leq c \leq f(b)$ gibt es ein $x_0 \in [a,b]$ so dass $f(x_0)=c$.\\
Für $c=0$ ist dieser Satz bei der Nullstellensuche hilfreich. \\
Suche Funktionswerte mit unterschiedlichem Vorzeichen: $f(a) \cdot f(b) < 0$.
Dann gibt es zwischen $a$ und $b$ mindestens eine Nullstelle.

\subsection{Bisektionsverfahren}
\begin{align*}
f(a) \cdot f(b)<0 \\
\hat{=} \text{Nullstellen in } (a,b) 
\end{align*}
Berechne Vorzeichen von $f(\frac{a+b}{2})$ \\
$\rightarrow$ $f(x)=0$ in $(a, \frac{a+b}{2})$ oder $(\frac{a+b}{2},b)$
$\rightarrow$ Berechne Vorzeichen von $f(\frac{a+b}{4})$ oder $\frac{3}{4}(a+b)$...\\
\\
\textbf{Beispiel:}\\ \\
$f(x)=x^3-x+0,3=0$\\
Wie viele Nullstellen?
\\
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline 
x & -2 & -1 & 0,5 & 1 \\ 
\hline 
f(x) & -5,7 & 0,3 & -0,075 & 0,3 \\ 
\hline 
\end{tabular} 
\end{center}
Wo sind die Nullstellen?\\
Bestimme die Nullstelle zwischen $x=0$ und $x=0,5$ genau auf eine Stelle nach dem Komma.
\begin{align*}
f(0,25)>0 \rightarrow \text{Nullstelle in } [0,25;0,5]\\
f(0,375)<0 \rightarrow \text{Nullstelle in } [0,25;0,375]\\
f(0,3125)>0 \rightarrow \text{Nullstelle in } [0,3125;0,375]
\end{align*}
Also ist die Nullstelle bei $0,3...$.

\subsection{Fixpunkt-Iteration}
Eine Gleichung der Form $x^{n+1}=f(x^{(n)})$ wird als Fixpunktgleichung bezeichnet.
Die Lösung(en) $\bar{x}$ mit $\bar{x}=f(\bar{x})$ heißen Fixpunkte. (da unter der Abbildung der Punkt $\bar{x}$ frei (=unveränderlich) bleibt.)
Jedes Nullstellenproblem kann als eine solche Fixpunktgleichung definiert werden.\\
\\
\textbf{Beispiel:}\\
Finde Nullstelle von $g(x)=x^3-x+0,3$.
Umformen der Fixpunktgleichung: $x^3-x+0,3=0$:
\begin{align*}
x^3-x+0,3=0 \rightarrow x=x^3+0,3 \\
\rightarrow x^{(n+1)}=(x^{(n)})^3+0,3
\end{align*}

Wir wählen zunächst Startwerte nahe der vermuteten Nullstellen ($\hat{=}$ Fixpunkten) und berechnen Werte für folgende n.\\
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline 
$n$ & $x^{(n)}$ & $x^{(n)}$ & $x^{(n)}$ \\ 
\hline 
0 & -1 & 0 & 1 \\ 

1 & -0,7 & 0,3 & 1,3 \\ 

2 & -0,043 & 0,327 & 2,497 \\ 

3 & 0,2999 & 0,3349 & 15,87 \\ 

4 & 0,327 & 0,337 &  \\ 

8 & 0,33877 & 0,33891 &  \\ 
\hline 
\end{tabular} 
\end{center}
Die ersten zwei Startwerte konvergieren zum selben Fixpunkt.
In der Tat ist nur der Fixpunkt $\bar{x}=0,3389...$ anziehend, die anderen werden als abstoßend bezeichnet.
Der Kreuzungspunkt zwischen der linken und der rechten Seite der Gleichung liefert den Fixpunkt.\\
Vergleich der Steigungen  von $y=x$ und $y=f(x)$ am Fixpunkt:\\
Steigung von $f(x)$ ist kleiner als $x$, also $f'(x)<1 \rightarrow$ Grund für Konvergenz.
Die Konvergenz ist also umso schneller je kleiner $f'(x)$ am Fixpunkt.\\
\\
\textbf{Fixpunktsatz:}\\ 
Sei $f:[a,b]\rightarrow \mathbb{R}$ mit stetiger Ableitung $f'(x)$ und $\bar{x}$ ein Fixpunkt von $f$, dann gilt für die Iteration:
\begin{align*}
x^{(n+1)}=(x^{(n)})^3+0,3:
\end{align*}
ist $|f'(\bar{x})|<1$ so konvergiert $x^{(n)}$ gegen $\bar{x}$ falls $x^{(0)}$ nage genug an $\bar{x}$ liegt: $\bar{x}$ ist ein anziehender Fixpunkt.\\
Ist $|f'(\bar{x})|>1$ so konvergiert $x^{(n)}$ für keinen Startwert $x^{(n)} nicht \bar{x} $. $\bar{x}$ ist dann ein abstoßender Fixpunkt.
\\
\\
\textit{zurück zum Beispiel:}\\
Plot der Funktionen und Kontrolle der Startpunkte.\\
Die Fixpunkte sind $\bar{x_1}=-1,125$, $\bar{x_2}=0,3389$ und $\bar{x_3}=0,7864$.\\
Plot der Abbildung $\rightarrow$ stabiler Punkt ablesbar.

\subsection{Newton Verfahren}
Gegeben: differenzierbare Funktion $f(x)$\\
Gesucht: Nullstelle $\bar{x}$ mit $f(\bar{x})=0$\\
Ausgangspunkt $x_0$ (in der Nähe von $\bar{x}$)\\
\textbf{Lösung:}\\
Linearisierung von $f(x)$ um $x_0$:\\
$f(x) \approx f(x_0)+(x-x_0) \ f'(x_0) = 0$\\
$f'(x_0) \neq 0$
$\rightarrow x_0 - \frac{f(x_0)}{f'(x_0)}$\\
Das heißt die Funktion $f(x)$ wird durch die Tangente am Pukt $x_0$ genähert. Verbesserungen sind im Prinzip möglich.\\
Wird dieses Prinzip iterativ angewandt redet man vom Newton-Verfahren.
\begin{align*}
x_{n+1} =x_{n} - \frac{f(x_n)}{f'(x_n)}
\end{align*}
An unserem Beispiel ergibt sich:
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline 
$n$ & $x_n$ &  &  \\ 
\hline
0 & -1 & 0 & 1 \\ 

1 & -1,15 & 0,3 & 0,85 \\ 

2 & -1,12615 & 0,33699 & 0,7951 \\ 

3 & -1,1254 & 0,3389 & 0,78668 \\ 

4 & -1,12542 & 0,33894 & 0,78649 \\ 
\hline 
\end{tabular} 
\end{center}
Nach nur 4 Iterationen hat man eine Genauigkeit von $10^{-4}$ erreicht.
Das Newton-Verfahren ist sehr schnell und beliebt; ein Nachteil liegt allerdings darin, dass in jedem Schritt eine Ableitung berechnet werden muss.\\
\\
\textbf{Lösung 1:} \textit{Das Vereinfachte Newton-Verfahren}\\
Statt in jedem Schritt $f'(x_n)$ zu berechnen wird immer wieder $f'(x_0)$ verwandt:
\begin{align*}
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_0)}
\end{align*}
Allerdings konvergiert dieser Ansatz nicht so schnell.\\
\\ 
\textbf{Lösung 2:} \textit{Sekantenverfahren} \\
Statt der Steigung im Punkt $x_n$ wird die Aleitung durch Differenzenbildung berechnet:
\begin{align*}
x_2=x_1-\frac{f(x_1) \ (x_1-x_0)}{f(x_1)-f(x_0)} \\
\rightarrow x_{n+1}=x_n-\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})} \cdot f(x_n)
\end{align*}
Hier wird also ebenfalls keine Ableitung benötigt; allerdings werden zwei Startwerte benötigt.
Die Konvergenzgeschwindigkeit ist nicht ganz so gut wie beim Newton-Verfahren, allerdings ist der Rechenaufwand gering. (In jedem Schritt muss nur eine Funktionsauswertung vorgenommen werden.)
Im Allgemeinen ist das Sekantenverfahren besser als das Newton-Verfahren.

\subsection{Konvergenzkriterien}
Die Effizienz von Nullstellensuche hängt von der Konvergenzgeschwindigkeit ab. \\
\textbf{Definition:} Sei $x_n$ eine Folge mit $lim_{n \rightarrow \infty}=\bar{x}$, dann hat die Folge eine Konvergenzordnung $q \geq 1$ wenn es die Konstante $c>0$ gibt mit:
\begin{align*}
|x_{n+1}| \leq c \ |x_n-\bar{x}|^q \text{für alle n}
\end{align*}
Für $q=1$ muss auch $c<1$ gelten.
Aus $q=1$ folgt die lineare Konvergenz aus $q=2$ die quadratische.\\
\\
\textbf{Beispiel:}\\
Sei $x_n$ eine lineare, konvergente Folge und $y_n$ eien quadratisch, konvergente Folge mit den Grenzwerten $x_n$ und $\bar{x}$.
Wir starten mit dem gleichen Abstand zur NS: $|x_n-\bar{x}|\leq 0,1$, $|y_n-\bar{x}|\leq 0,1$.
Im nächsten Schritt: $|x_{n+1}-\bar{x}| \leq c \cdot 0,1$ und $|y_{n+1}-\bar{x}| \leq c \cdot 0,01$.
Das quadratische Verfahren konvergiert also deutlich schneller.\\
Bemerkung: Für einfache Nullstellen konvergiert das Newtonverfahren quadratisch, das Sekantenverfahren mit $q=\frac{\sqrt{5}+1}{2}$ und das vereinfachte Newton-Verfahren linear.

Außerdem lässt sich für mehrfache Nullstellen das Newton-Verfahren verbessert werden und eine quadratische Konvergenz erreicht werden:
\begin{align*}
x_{n+1}=x_n + n \ \frac{f(x_n)}{f'(x_n)},
\end{align*}
wobei $m$ die Vielfachheit der Nullstelle beschreibt.

Cave: numerische Auslöschung bei $f'(x_n) \approx 0$.

\subsection{Zusammenfassung}
\begin{itemize}
\item[Bisektion:] Einfaches Verfahren, schlechte Konvergenz, Fehler halbiert sich mit jedem Iterationsschritt.
Allerdings wird dieses Verfahren häufig verwendet um die Nullstelle einzugrenzen und dann ein Verfahren mit besserer Konvergenz zu nutzen.
\item[Fixpunktiteration:]Linearkonvergenz bei anziehenden Nullstellen, aber Vorsicht bei abstoßenden Nullstellen.
Diese Methode ist allerdings einfach anzuwenden und numerisch günstig, da immer nur eine Rechenoperation ausgeführt werden muss.
\item[Newton-Verfahren:] Quadratische Konvergenz bei einfachen Nullstellen, allerdings die die Berechnung aufwendig, da ABleitungen gebildet werden müssen.
Ein Vorteil ist die Möglichkeit der Anwendung auf mehrdimensionale Probleme, allerdings benötigt dann jede Iteration die Lösung eines linearen Gleichungssystems.
\item[Vereinfachtes Newton-Verfahren:]Weniger Rechenaufwand, da die Ableitung nicht neu berechnet werden muss allerdings nur lineare Konvergenz.
\item[Sekantenverfahren:] Eine relativ gute Konvergenz mit $q\approx 1,618$ und auch gute numerische Effizienz, da nicht die Ableittung sondern die Differenzenquotient genutzt wird. Allerdings muss die Auslöschung bei Nullstellen beachtet werden: $\frac{f(x_n)}{f'(x_n)} \approx \frac{0}{0}$.
\end{itemize}

\section{Lineare Gleichungssysteme}
WIr betrachten ein Gleichungssystem $a_{11} x_1+a_{12} x_2+...+a_{1n} x_n=b_1$, $a_{21} x_1+a_{22} x_2+...+a_{2n} x_n=b_2$ und $a_{31} x_1+a_{32} x_2+...+a_{3n} x_n=b_3$. Ein solches Gleichungssytem ist exakt lösbar wenn die Anzahl der unabhängigen Gleichungen gleich der Anzahl der Unbekannten ist. 
Die Lösung von größeren Gleichungssystem erfordert numerische Verfahren, dabei unterscheidet man zwischen zwei verschiedenen Klassen.
\begin{itemize}
\item[Direkte Verfahren:] In einer endlichen Anzahl von Schritten erhält man die exakte Lösung im Rahmen der numerischen Genauigkeit.
\item[Iterative Verfahren:] Hier wird eine Folge $b_k$ von Lösungsvektoren erzeugt die gegen b konvergiert.
\end{itemize}

\subsection{Gauß-Verfahren}
Ein Gleichungssystem der Form $\begin{pmatrix}
a_{11} & a_{12} & ... & a_{1n} \\
0 & a_{22} & ... & a_{2n} \\
. & .\\
. & .\\
. & .\\
0 & 0 & ... & a_{nn}
\end{pmatrix}	\begin{pmatrix}
x_1 \\ x_2 \\. \\ . \\ . \\x_n  
\end{pmatrix}=\begin{pmatrix}
b_1 \\ b_2 \\. \\. \\. \ \\x_n  
\end{pmatrix}$

Um eine beliebige Matrix auf diese dieser Form zu bringen können die Zeilen einzeln durch elementare Zeilenoperationen umgeformt werden:
\begin{align*}
z_j=z_j-\frac{a_{ji}}{a_{ii}},
\end{align*}
für Spalte $i$ und Zeile $j$ mit $i=1,..,n$ und $j=i+1,...,n$.

Das Ergebnis lässt sich dann durch Rückeinsetzung lösen.\\
Bemerkung: DIe Anzahl der benötigten Schritte wächst beim Gauß-Verfahren kubisch mit der Größe des LGS $n$: $\sigma(n^3)$.\\
\\
\textbf{Beispiel: Fehlerminimierung mit Pivotisierung}
$A= \begin{pmatrix}
-16^{-4} & 1 \\
2 & 1 
\end{pmatrix}$ und $b=\begin{pmatrix}
1 \\ 0
\end{pmatrix}$, wir rechnen in 4-stelliger Dezimalgenauigkeit.
Aus dem Gauß-Algorithmus folgt dann:
\begin{align*}
z_2=[0,20000+1,20000]
\end{align*}
Beim Rückeinsetzen ergibt sich dann:
\begin{align*}
x_2=\frac{20000}{20000+1} \approx 1 \\
x_1=\frac{1}{-10^{-4}} \cdot (1-1) =0.
\end{align*}
Die exakten Lösungen wären allerdings $x_1=-0,499975$ und $x_2=0,99995$.

Um einen kleineren Fehler zu erhalten werden zunächst die erste und die zweite Zeile getauscht. Daraus ergibt sich dann:
\begin{align*}
z_2=[-10^{-4+10^{-5}} \approx 0,1+0,5\ 10^{-4} \approx 1, 1+0]\\
x_2=1 \\
x_1=- \frac{1}{2}.
\end{align*}
\textbf{Satz:} Der Fehler im Gauß-Verfahren wird durch die sogenannte \textit{Spaltenpivotisierung} minimiert. Dabei werden vor jedem Eliminationsschritt für die i-te Spalte die Zeilen des LGS so umsortiert, dass gilt:
\begin{align*}
|a_{ii}|= \text{max}\{|a_{ij}|, j=1,...n \}.
\end{align*}
Dann gilt im Eliminationsschritt für den Fehler. $|\frac{a_{ji}}{a_{ii}}|\leq 1$.
Dann wird der Gauß-Algorithmus für die Spalte i angewandt und für $i+1$ erneut auf das betragsmäßig größte Element in den Zeilen $j=i+1,...n$
geprüft.

\subsection{LR(LU)-Zerlegung}
Möchte man mehrere LGS mit der selben Koeffizientenmatrix aber anderer rechten Seite lösen so empfiehlt es sich die Elementaren Umformungen zu merken. \\ 
Für $A= \begin{pmatrix}
1  & 2 & 3\\
6  & -2 & 2\\
-3  &1 & 4
\end{pmatrix}$
führten wir durch:\\
\begin{align*} 
 z_2=z_2-6z_1 \rightarrow L_1=\begin{pmatrix}
1  & 0 & 0\\
-6  & 1 & 0\\
0  &0 & 1
\end{pmatrix} \\
z_3=z_3+3z_1 \rightarrow L_2=\begin{pmatrix}
1  & 0 & 0\\
0  & 1 & 0\\
3  &0 & 1 
\end{pmatrix} \\
z_3=z_3+\frac{1}{2} z_2 \rightarrow L_1=\begin{pmatrix}
1  & 0 & 0\\
0  & 1 & 0\\
0  &\frac{1}{2} & 1
\end{pmatrix}.
\end{align*}

Daraus folgt dann:
\begin{align*}
L=L_3 \ (L_2 \cdot L_1)=\begin{pmatrix}
1  & 0 & 0\\
-6  & 1 & 0\\
0  &\frac{1}{2} & 1
\end{pmatrix}.
\end{align*}
Damit ergibt sich dann die rechts-obere Dreiecksmatrix R zu:
\begin{align*}
R=L(A)
\end{align*}
So lässt sich dann für ein beliebiges LGS $A \ x = c$ der Lösungsvektor aus $L(c)$ berechne denn da die linke Seite von
\begin{align*}
R=L(c)
\end{align*}
bereits eine rechts-obere Dreiecksmatrix ist kann die Gleichung direkt durch Rückeinsetzen gelöst werden. \\
\\
\textbf{Satz} Für jede $nxn$-Matrix, für die der Gauß-Algorithmus ohne Zeilenvertauschungen durchführbar ist, gibt es $nxn$-Matrizen $L$ und $R$ mit den Eigenschaften:
\begin{itemize}
\item $L$ ist eine links-untere Dreiecksmatrix mit $l_{22}$ für $i=1,...,n$
\item R ist eine rechts-obere Dreiecksmatrix mit $r_{22} \neq 0$ für $i=1,...,n$
\item $A=L^{-1} \ R$ bezeichnet man als LR-Zerlegung von $A$. 
\end{itemize}
Es gilt:
\begin{align*}
A \ x = b \Leftrightarrow L b = y \quad \text{und} \quad R x =y
\end{align*}
\subsection{Cholesky-Zerlegung}
\textbf{Definition:} Eine symmetrische $n x n$-MatrixA heißt positiv-definit wenn gilt:
\begin{align*}
x^T Ax>0
\end{align*}
\textbf{Satz}Für jede positiv-definite Matrix $A$ gibt es genau eine rechts-obere Dreiecksform mit $r_{ii}>0 \quad \text{für} \quad 1,...,n$
und $A=R^TR$, diese Zerlegung nennt man Cholesky-Zerlegung. 

Die Berechnung der Zerlegung erfolgt wie folgt:\footnote{http://de.wikipedia.org/wiki/Cholesky-Zerlegung\#Pseudocode} 
\begin{lstlisting}[frame=single]


    For i = 1 To n
        For j = 1 To i
            Summe = a(i, j)
            For k = 1 To j-1
                Summe = Summe - a(i, k) * a(j, k)
            If i > j Then
                a(i, j) = Summe / a(j, j)   
            Else If Summe > 0 Then          
                a(i, i) = Sqrt(Summe)       
            Else
                ERROR                       
\end{lstlisting}

Damit ergibt sich für $A=\begin{pmatrix}
4 & 4 & 2 \\
4 & 5 & 5 \\
4 & 5 & 26 \\
\end{pmatrix}$, $R=\begin{pmatrix}
2 & 2 & 1 \\
0 & 1 & 3 \\
0 & 0 & 4 \\
\end{pmatrix}$.\\
\textbf{Bermerkung} Der numerische Aufwand des Verfahrens beträgt $(\frac{1}{6} n^3 + \frac{1}{2} n^2 -\frac{2}{3} n)$ und eine Wurzelberechnungen.
Das Gauß-Verfahren benötigt$(\frac{n^3}{3}+\frac{n}{3})$, ist also für $n \geq 2 $ langsamer.

\subsubsection{Fehlerrechnung bei LGSW}
WIr wollen untersuchen wie sich Fehler in einem LGS auf dessen Lösung auswirken, dazu benötigen wir ein Maß.
\textbf{Definition:} Eine Abbildung $||\ ||: \mathbb{R}^n \rightarrow \mathbb{R}$ heißt Vektornorm, wenn für alle $x,y \in \mathbb{R}^n$ und alle $\lambda \in \mathbb{R}$ gilt:

\begin{itemize}
\item $||x|| \geq 0$ und $||x||=0$ wenn $x=0$
\item $|| \lambda x|| = |\lambda| ||x||$
\item $||x+y|| \leq ||x||+||y||$
\end{itemize}

Beispiele hierfür sind Summennorm, Maximunmnorm und die euklidische Norm.
Diese sind äquivalent wenn gilt:...\\

Entsprechend lassen sich auch Matrixnormen definieren: Spaltensummennorm Spektralnorm, Zeilensummennorm.

\textbf{Bemerkung:} Matrixnormen erfüllen die Eigenschaften der zugrundliegenden Vektornormen, insbesondere auch die Äquivalenz:
\begin{align*}
||A x||_v \leq ||A||_v||x||_v,
\end{align*}
man sagt dann: die Norm ist kompatibel.\\
\textbf{Satz} Sei $|| \ ||$ eine Norm und A eine reguläre $nxn$ Matrix und:$A \ x=b$ und $A \ x' =b'$

Dann gilt:
\begin{align*}
A(x-x')=b-b' \\
x-x'=A^{-1}(b-b')\\
||x-x'|| \leq ||A^{-1}|| \ |b-b'|.
\end{align*}
Und ensprechend mit der MUltiplikation von $||A|| ||x|| \geq ||b|| \rightarrow ||\frac{1}{x}|| \leq \frac{||A||}{||b||}$:
\begin{align*}
\frac{||x-x'||}{||x||} \leq ||A|| ||A^{-1}|| \cdot \frac{||b-b'||}{||b||}.
\end{align*}
Man nennt $||A|| \ ||A^{-1}||=cond(A)$ die Konditionszahl der Matrix $A$ bzgl. der verwendeten Norm. \\
\textbf{Bermerkung:} Die Koordinationszahl gibt also die max. Verstärkung des relativen Fehleran, währen $||A^{-1}||$ die maximale Verstärkung des absoluten Faktors ist.\\
\textbf{Satz:} Ist nicht nur die rechte Seite eines LGS fehlerbehaftet, sondern auch die Koeffizientenmatrix, so gilt für die Lösungen der beiden LGS $A \ x = b$ und $A' \ x' = b'$ mit $\Delta A= A-A'$ und $\Delta x= x-x'$
\begin{align*}
\frac{||\Delta x||}{||x||} \leq \frac{cond(A)}{1-cond(A) \frac{||\Delta|A|}{A}}   (norm( deltA(A + norm( deltab/b)
\end{align*}

\subsection{Iterative Verfahren}
\textbf{Idee:} Das Linerares Gleichungssystem $A x=b$ ist äquivalent zu einem vektoriellen Nullstellenproblem:
\begin{align*}
Ax-b=0.
\end{align*}
Anstelle des Nullstellenproblems betrachtet man dann das Fixpunktproblem.
Um auf der rechten Seite der Matrix ein $x$ zu erzeugen muss die Matrix zerlegt werden: $A=I+A-I$, mit der Einheitsmatrix $I$. Daraus folgt:
\begin{align*}
0=(I+A-I)x-b \Rightarrow I x=(I-A)x+b\\
\Rightarrow x=(I-A)x+b.
\end{align*}
Dies entspricht dann einer vektoriellen Fixpunktgleichung, welche als Iterationsgleichung aufgefasst werden kann.

Mit einem Startwer $x^{(0)}$ folgt damit:
\begin{align*}
x^{(1)}=(I-A)x^{(0)}+b, \\
x^{(n+1)}=(I-A)x^{(n)}+b.
\end{align*}
Ein Fixpunkt der Itertaionsgleichung einspricht dann also der Lösung des LGS.

\textbf{Beispiel:}\\
Aus $A=\begin{pmatrix}
4 & -1 & 1 \\
-2 & 5 & 1 \\
1 & -2 & 5 \\
\end{pmatrix}$, 
$b=\begin{pmatrix}
5\\
11 \\
12\\
\end{pmatrix}$ folgt:
$x^{(n+1)}=\begin{pmatrix}
-3 & 1 & -1 \\
2 & -4 & -1 \\
-1 & 2 & -4 \\
\end{pmatrix} x^{(n)}+
\begin{pmatrix}
5\\
11 \\
12\\
\end{pmatrix} $,
was nicht konvergiert.

\subsection{Jacobi Verfahren(Gesamtschnitt)}
Betrachte Zerlegung $A=\underbrace{\begin{pmatrix}
0 & 0 & 0 & .\\
a_{21} & 0 & 0 & . \\
a_{31} & a_{32} & 0 & .\\
.&.&.&.
\end{pmatrix}}_L+ \underbrace{\begin{pmatrix}
a_{11} & 0 & 0 & .\\
0 & a_{22} & 0 & . \\
0 & 0 & a_{33} & .\\
.&.&.&.
\end{pmatrix}}_D+
\underbrace{\begin{pmatrix}
0 & a_{12} & a_{13} & .\\
0 & 0 & a_{23} & . \\
0& 0 & 0 & .\\
.&.&.&.
\end{pmatrix}}_R$.
Daraus ergibt sich dann:
\begin{align*}
Ax-b=0=(L+D+R)x-b\\
\Rightarrow Dx=-(L+R)x+b \\
\Rightarrow x=-D^{-1}(L+R)x+D^{-1}b \\
\Rightarrow x^{(n+1)}=-D^{-1}(L+R)x^{(n)}+D^{(-1)}b
\end{align*}
Dies liefert für unser Beispiel von vorher:
\begin{align*}
x^{(n+1)}=\begin{pmatrix}
0 & 0,25 & -0,25 \\
0,4 & 0 & -0,2 \\
-0,2 & 0,4 & 0 \\
\end{pmatrix}x^{(n)}+\begin{pmatrix}
1,25\\
2,2\\
2,4
\end{pmatrix}.
\end{align*}
Dies liefert dann:
\begin{tabular}{c|cccccc}

$i$ & 0 & 1 & 2 & 3 & 4 & 5 \\ 
\hline 
$x^{(i)}$ & 0 & 1,25 & 1,2 & 1,0475 & 1,0065 & 0,9973 \\ 
 
& 0 & 2,2 & 2,22 & 2,074 & 2,0094 & 1,9986 \\ 

& 0 & 2,4 & 3,03 & 3,048 & 3,0201 & 3,0024 \\ 

\end{tabular} 

dies konvergiert gegen: $\bar{x}=\begin{pmatrix}
1\\
2\\
3
\end{pmatrix}$.

\subsubsection{Einzelschrittverfahren (Gauß-Seidel Verfahren)}
Jacobi-Verfahren in Komponenten:
\begin{align*}
x_1^{(n+1)}=0,25 x^{(n)}_2-0,25 x^{(n)}_3 +1,25\\
x_2^{(n+1)}=0,4 x^{(n)}_1-0,2 x^{(n)}_3 +2,2\\
x_3^{(n+1)}=-0,2 x^{(n)}_1+0,4 x^{(n)}_2 +2,4.
\end{align*}
Wenn man annimmt, dass $x^{(n+1)}$ komponentenweise näher am Lösungsvekotr liegt, sollte man in der zweiten Gleichung $x_1^{n+1}$ u.s.w. benutzen, damit ergibt sich:
\begin{align*}
x_1^{(n+1)}=0,25 x^{(n)}_2-0,25 x^{(n)}_3 +1,25\\
x_2^{(n+1)}=0,4 x^{(n+1)}_1-0,2 x^{(n)}_3 +2,2\\
x_3^{(n+1)}=-0,2 x^{(n+1)}_1+0,4 x^{(n+1)}_2 +2,4.
\end{align*}
In Matrix Schreibweise:
\begin{align*}
x^{(n+1)}=\begin{pmatrix}
0 & 0,25 & -0,25\\
0& 0 & -0,2\\
0& 0 & 0
\end{pmatrix} x^{(n)}+
\begin{pmatrix}
0 & 0 & 0\\
0,4& 0 & -0\\
-0,2& 0,4 & 0
\end{pmatrix} x^{n+1}+
\begin{pmatrix}
1,25\\
2,2\\
2,4
\end{pmatrix},
\end{align*}
oder einfacher:
\begin{align*}
x^{(n+1)}=-D^{-1}(R x^{(n)}+L x^{(n+1)}-b)\\
\Leftrightarrow (D+L)x^{(n+1)}=-R x^{(n)}+b\\
\Leftrightarrow x^{(n+1)}=-(D+L)^{-1} R x^{(n)}+(D+L)^{-1} b.
\end{align*}
In der Praxis wird nicht das Inverse der Matrix $(D+L)$ berechnet, sondern das Gleichungssystem $(D+L)x^{(n+1)}=-R x^{(n)}+b$ wird gelöst(durch Vorwärtseinsetzen).\\
Zurück zum Beispiel:
\begin{tabular}{c|ccccc} 
 & 0 & 1 & 2 & 3 & 4 \\ 
\hline 
 & 0 & 1,25 & 1,1175 & 1,006 & 1,001 \\ 
 
 & 0 & 2,7 & 2,001 & 2,007 & 2,0002 \\ 

 & 0 & 3,24 & 2,977 & 3,002 & 2,9998 \\ 

\end{tabular},
konvergiert also schneller.\\
\subsubsection*{Konvergenz von linearen Iterationsverfahren} 
Gegeben sei eine lineare Matrix-Fixpunkt-Iterationsgleichung $x^{n+1}=B x^{(n)}+b$, wobei B eine $m x m$ Matrix, $b \in \mathbb{R}^m$ und $\bar{x}$ ein Fixpunkt mit: $\bar{x}= B \bar{x}+ b$.\\
Sei $|| . ||$ eine Matrixnorm, dann gilt: \\
$\bar{x}$ anziehender Fixpunkt falls $||B|| < 1$\\
$\bar{x}$ abstoßender Fixpunkt falls $||B|| > 1$\\
\\
\textbf{Matrix-Fixpunktsatz für anziehenden Fixpunkt:} Fixpunktiteration konvergiert für alle Startwerte.\\
\textit{Bemerkungen:} F+r Gesamtschrittverfahren (Jacobi) gilt \textit{$B=-D^{-1}(L+R)$}. Für die $\infty$-Norm (Zeilensummennorm)gilt:
\begin{align*}
||B||_{\infty}=max_{i,m} \sum_{j \neq i} \frac{|a_{ij}|}{|a_{ii}|}=max_{i,m} \frac{1}{|a_{ii}|} \sum_{j \neq i} |a_{ij}|<1\\
\Leftrightarrow |a_{ii}|>\sum_{j \neq i} |a_{ij}|,
\end{align*}
für alle $i$, dies wird Zeilensummenkriterium genannt.
Eine solche Matrix heißt diagonaldominant.

Im Einzelschrittverfahren ist $B=-(D+L)^{-1} R$ und man kann zeigen, dass $||(D+L)^{-1}R||_{\infty} \leq ||D^{-1}(L+R)||_{\infty}$.
Konvergiert also das Gesamtschrittverfahren so konvergiert auch dasd Einzelschrittverfahren, die Umkehrung gilt nicht.

\subsection{Nichtlineare Gleichungssyteme-Newton Verfahren}
häufig: $n$ nicht-lineare Gleichungen mit $n$ Unbekannten\\
gegeben: vektorielle Funktion  $f:\mathbb{R}^n \rightarrow \mathbb{R}^n$\\
gesucht: Vektor $\bar{x} \in \mathbb{R}^n$ mit $f(\bar{x})=0$.\\
$f(x)=f(x_1,...,x_n)= \begin{pmatrix}
f_1(x_1,...,x_n) \\
f_2(x_1,...,x_n) \\
...\\
f_3(x_1,...,x_n)
\end{pmatrix}= \begin{pmatrix}
0\\
0\\
...\\
0
\end{pmatrix}
$.
Es gibt kein einfaches Verfahren um zu prüfen ob ein Gleichungs-System lösbar ist und wieviele Lösungen es gibt.
$\rightarrow$ verallgemeinere Newton-Verfahren auf $n$-Dimensionen.
\begin{align*}
g: \mathbb{R} \rightarrow \mathbb{R} \Rightarrow g(x) \approx g(x_0)+ (x-x_0)g'(x_0)=0\\
f: \mathbb{R}^n \rightarrow \mathbb{R}^n \Rightarrow f(x) \approx g(x^{(0)})+ (x-x^{(0)})D(f(x^{(0)}))=0
\end{align*}
Zu lösen ist also das Gleichungssystem:
\begin{align*}
x=x^{(0)}-(D(f(x^{(0)})))^{-1} f(x^{(0)}) \\
x^{(n+1)}=x^{(n)}-(D(f(x^{(n)})))^{-1} f(x^{(n)}).
\end{align*}
Dies wird als Newton-Verfahren bezeichnet.

\subsubsection*{Jacobi Matrix}
Die Jacobi Matrix der partiellen Ableitungen ist definiert als:
\begin{align*}
D(f(x))=\frac{\partial f_i(x)}{\partial x_j}=\begin{pmatrix}
\frac{\partial f_1(x)}{\partial x_1} & \frac{\partial f_1(x)}{\partial x_2} & ... & \frac{\partial f_1(x)}{\partial x_n} \\
.\\
.\\
.\\
\frac{\partial f_n(x)}{\partial x_1} & ... & & \frac{\partial f_n(x)}{\partial x_n} 
\end{pmatrix}.
\end{align*}
So ergibt sich mit $f(x_1,x_2)= \begin{pmatrix}
2 x_1 + 4 x_2 \\ 4 x_1++ 8x_2^3
\end{pmatrix}$ und $x^{(0)}= \begin{pmatrix}
4 \\ 2
\end{pmatrix}$ $\bar{x}= \begin{pmatrix}
-2 \\
1
\end{pmatrix}$

\textbf{Satz:} Das Newton Verfahren konvergiert quadratisch wenn:
\begin{itemize}
\item $x^{(0)}$ nahe genug an $\bar{x}$ liegt
\item $D(f(\bar{x}))$ regulär ist
\item $f$ dreimal stetig differenzierbar ist
\end{itemize}

\textbf{Achtung:} Nicht immer konvergiert das Newton Verfahren gegen eine Nullstelle. Ein Beispiel hierfür ist die Funktion: $f= \begin{pmatrix}
x_1^3 -x_2 -1 \\
x_1^2-x_2
\end{pmatrix}$.
Die liegt daran, dass die zugehörige Jacobi-Matrix nicht regulär (nicht invertierbar) ist.

\textbf{Bemerkung:} in der Praxis werden partielle Ableitungen durch Differenzieren angenähert: $D(f(x))_{i j}=\frac{f_i(x+h \ e_j)- f_i(x)}{h}$,
mit dem Einheitsvektor in $j$-Richtung $e_j$.


\subsubsection*{Vereinfachtes Newton-Verfahren}
Der Rechenaufwand pro Schritt wird dadurch minimiert, dass stets $D(f(x^{(0)}))$ benutzt wird, dadurch ist die Konvergenz nur noch linear.

\section{Interpolation und Ausgleichsrechnung}
In vielen Anwendungen: Messdaten sollen durch eine Formel beschrieben werden, zum Beispiel um Integrale, Differentiale, etc. zu berechnen.\\
\textbf{Definition:} gegeben seien $n+1$ Wertepaare $(x_1, f_i)$ mit $i=(0,...,n)$; gesuchten stetigen Funktionen $f(x)$ mit $f(x_i)=f_i$ für alle $i=0,..,n$.\\
$\{x_n\}$: Stützstellen, $\{f_n\}$: Stützwerte; $\{x_n, f_n\}$ Stützpunkte und $f(x)$ Interpolierende der Stützpunkte.\\
\textbf{Bemerkung:} Das Interpolationsproblem ist nicht eindeutig.

\subsubsection{Polynom-Interpolation}
Ein Polynom $n$-ten Grades besitzt $n+1$ Freiheitsgrade, es ist also möglich die $n+1$ Koeffizienten des Polynoms so zu wählen, dass $f(x_i)=f_i$ ist!\\
\textbf{Satz:} Gegeben seien $n+1$ Wertepaare $(x_i,f_i)$, dann gibt es genau ein Polynom vom Grad höchstens $n$ so dass das Polynom $p(x_i)=f_i$ für alle $i$.
Hierbei spricht man bei $n=1$ von einer linearen, bei $n=2$ von einer quadratischen Interpolation und so weiter.\\
\subsubsection*{ Ineffiziente Methode:}
Explizite Lösung eines linearen Gleichungssystems.\\
\textbf{Beispiel:}\\
\begin{figure}
\center

\begin{tabular}{|c||c|c|c|c|}

$x_i$ & -1 & 0 & 1 & 2 \\ 
\hline 
$f_i$ & 5 & -2 & 9 & -4 \\ 

\end{tabular}\\  
\end{figure}
\begin{align*}
p(x)=\sum_{i=0}^3 a_i x^i \\
f_j=\sum_{i=0}^3 a_i x_j^i,
\end{align*}
also 4 Gleichungen für 4 Unbekannte:
\begin{align*}
5=a_0+a_i(-1)+a_2 (-1)^2 + a_3 (-1)^3 \\
-2=a_0 \\
9=a_0+a_1 (1) + a_2 (1)^2 + a_3 (1)^3\\
-4 =a_0 + a_1 (2) + a_2 (2)^2 +a_3 (2)^3.
\end{align*}
Als Lösung ergibt sich: \begin{align*}
a_0=2, \quad a_1=9, \quad a_2=9 \\
\Rightarrow p(x)=-2+9 x + 9 x^2 - 7 x^3.
\end{align*}

\subsubsection*{Lagrange Form}
\begin{align*}
p(x)=\sum_{i=0}^3 f_i \ l_i(x)
\end{align*}
mit $l_i(x)=\prod_{j=0, j \neq i}^n \frac{x-x_j}{x_i-x_j}$. \\
In unserem Beispiel: $l_0(x)=(\frac{x-x_1}{x_0-x_1})(\frac{x-x_2}{x_0-x_2})(\frac{x-x_3}{x_0-x_3})$. Hier ist zu sehen, dass $l_0(x)$ an allen Stützstellen außer an $x_0$, und $l_0(x_0)=1$. (analoges gilt für alle anderen Polynome $l_1$, $l_2$, $l_3$)
$\rightarrow p(x_i)=f_i$. \\
Jedes $l_1$ ist ein Polynom dritten Grades.
Die Berechnung von $p(x)$ durch Summation von $l_i(x)$ ist allerdings immer noch aufwendig, es gibt eine schnellere Methode.

\subsubsection*{Newtonsche Interpolationsformel}

Gegeben sind Wertepaare $(x_i,f_i)$ für $i=0,...n$.
Zunächst werden "dividierte Differenzen" berechnet:\\
(Pseudocode hier)\\
\textbf{Bemerkung:} $x_i$ müssen nicht nach Größe sortiert sein außerdem können neue Wertepaare angefügt werden.

In unserem Beispiel ergeben sich die dividierten Differenzen für $k=1$ zu:
\begin{align*}
f(x_0,x_1)=-7\\
f(x_1,x_2)=11\\
f(x_2,x_3)=-13
\end{align*}
und für für $k=2$:
\begin{align*}
f(x_0,x_1,x_2)=9 \\
f(x_1,x_2,x_3)= -12
\end{align*}
Und daraus ergibt sich:
\begin{align*}
f(x_0,x_1,x_2,x_3)=-7.
\end{align*}

\subsubsection*{Interpolationsformel}
\begin{equation*}
p(x)= \sum_{i=0}^n f(x_0,...,x_i) \prod_{j=0}^{i-1} (x-x_j)
\end{equation*}
\begin{equation*}
=f(x_0)+f(x_0,x_1)(x-x_0)+f(x_0,x_1,x_2)(x-x_0)(x-x_i)
\end{equation*}
\begin{equation*}
+f(x_0,x_1,x_2,x_3)(x-x_0)(x-x_1)(x-x_2)
\end{equation*}

Für unser Beispiel ergibt sich also:
\begin{equation*}
p(x)=5+(-7)(x-(-1))+9(x-(-1))(x-0)+(-7)(x+1) x (x-1)
\end{equation*}
\begin{equation*}
=-7x^3+9x^2+9x-2
\end{equation*}

Außerdem kann p(x) durch Ausklammern umgeformt werden:
\begin{equation*}
p(x)=f(x_0)+(x-x_0)(f(x_0,x_1)+(x-x_1)(f(x_0,x_1,x_2)+(x-x_2) f(x_0,x_1,x_2,x_3))))
\end{equation*}
Durch das einsparen von Punktoperationen kann so Rechenzeit gespart werden.
\subsubsection*{Iterative Auswertung des Newton Polynoms, Horner Schema}
\begin{lstlisting}[frame=single]
	r_n=f(x_0,x_1,...,x_n)
        For k = n-1,0
       		r_k=r_(k+1) (x-x_k)+f(x_0,...,x_k)
       	End for p(x)=r_0                      
\end{lstlisting}
Man benötigt also zur Auswertung des Newton Polynoms nur n Punkt-Operatoren.
Daraus ergibt sich für unser Beispiel für:
\begin{align*}
r_3=f(x_0,x_1,x_2,x_3)=-7 \\
r_2=r_3 (x-x_2)+ f(x_0,x_1, x_2)=-7(x-1)+9 \\
r_1=r_2(x-x_1)+f(x_0,x_1)=(-7(x-1)+9)x-7 \\
r_0=r_1(x-x_0)+f(x_0)=((-7(x-1)+9)x-7)(x+1)+5
\end{align*}

\subsubsection*{Problem bei Polynom-Interpolation}
Gegeben sei die Funktion $f(x)=\frac{1}{1+x^2}$ mit: \\
\begin{tabular}{c|ccccccc}
 
$x_i$ & -3 & -2 & -1 & 0 & 1 & 2 & 3 \\ 
\hline 
$f_i$ & 0,1 & 0,2 & 0,5 & 1 &  &  &  \\ 

\end{tabular}
$\Rightarrow p(x)=1-0,64x^2+0,15x^4-0,06x^6$
Das Polynom weißt an den Rändern des Stützstellenintervalls Überschwinger auf.
Erhöhung von n verbessert die Situation nicht.
Außerdem divergiert das Polynom außerhalb des Intervalls. \\
$\rightarrow$ Interpolationspolynome sind nicht für Extrapolation geeignet.

\subsection{Spline Interpolation}
Ausgangsüberlegung: eine gute Näherung einer Funktion hat viele Stützstellen, allerdings besitzt ein Polynom hoher Ordnung auch viele Überschwingungen und divergiert stark.\\
\textbf{Idee:}Teile Intervall in Teilintervalle auf, und benutze jeweils Polynome niedriger Ordnung.\\
\begin{itemize}
\item \textbf{Einfachste Stufe:} Stückweise Interpolation mit Geraden, hierbei entstehen allerdings Knicke.\\
\item \textbf{Kubische Splines:} Hier wird zur ABhilfe eine stetige erste und zweite Ableitung gefordert.
\end{itemize}

\textbf{Definition:} eine Funktion $s(x):[a,b] \rightarrow \mathbb{R}$ heißt kubischer Spline zu den Stützstellen $a=x_0<x_1,...<x_n=b$ falls gilt:
\begin{itemize}
\item $s''(x)$ existiert und ist stetig
\item $s(x)$ ist auf den Intervallen $[x_i,x_{i+1}]$ jeweils ein Polynom 3. Grades
\item interpolieren der Spline Werten $f_0...f_n:s(x_i)=f_i \quad \text{für} \quad i=0,...,n$
\item periodischer Spline: zusätzlich $s(a)=s(b)$ und $s'(a)=s'(b)$ und $s''(a)=s''(b)$
\item natürlicher Spline: zusätzlich $s''(a)=s''(b)=0$
\end{itemize}
\subsubsection*{Bestimmung der Koeffizienten}
Ein Polynom 3. Grades hat 4 Koeffizienten $\hat{=}$ Freiheitsgrade; für n Intervalle haben wir also $4n$ Koeffizienten!

\begin{itemize}
\item Pro Intervall 2 Interpolationsbedingungen (rechts und links) $\rightarrow$ $2n$
\item An den Stützstellen $x_1,...,x_{n-1}$ sollen 1. und 2. Ableitungen gleich sein $\rightarrow$ $2(n-1)$ Bedingungen, zusammen gibt es also $4n-2$ Bedingungen.
\end{itemize}
Es sind also 2 zusätzliche Randbedingungen nötig.

\subsubsection*{Praktische Berechnung von kubischen Splines}
Wir definieren die "Momente" eines interpolierenden Splines als $M_i=s''(x_i)$. Mit diesen Momenten erhalten wir im Intervall $[x_{i-1},x_i], i=1,...,n$.
\begin{align*}
S_i(x)=M_{i-1} \frac{(x_i-x)^3}{6 h_i} + M_i \frac{(x-x_{i-1})^3}{6 h_i}+C_i(x- \frac{x_{i-1}+x_i}{2})+D_i\\
\rightarrow s'_i(x)= -M_{i-1} \frac{(x_i-x)^2}{2 h_i} + M_i \frac{(x-x_{i-1})^2}{2 h_i + C_i} \\
\rightarrow s''_i(x)= M_{i-1} \frac{x_i-x}{h_i} + M_i \frac{(x-x_{i-1})}{h_i}\\
h_i=x_i-x_{i-1}
\end{align*}
Die Ableitungen stimmen dabei an den Stützstellen überein:\\
z.B.: $s''_i(x=x_i)=s''_{i+1}(x=x_i) \Leftrightarrow M_i=M_i$ \\
$s'_i(x_i)=s'_{i+1}(x_i)$ und $s_i(x_i)=s_{i+1}(x_i)$ führen zu:
\begin{align*}
\mu_i M_{i-1} + 2 M_i + (1- \mu_i) M_{i+1} = 6 f(x_{i-1},x_i,x_{i+1})\\
C_i = \frac{f_i-f_{i-1}}{h i}- \frac{h_i}{6}(M_i-M_{i-1})\\
D_i= \frac{f_i+f_{i-1}}{2} -\frac{h_i^2}{12} (M_i+M_{i-1})\\
\mu_i = \frac{x_i-x_{i-1}}{x_{i+1}-x_{i-1}}
\end{align*}
Zusätzlich muss gelten:
\begin{itemize}
\item natürliche splines: $M_0=M_n=0$
\item preiodische splines: $M_0=M_n$, $\mu_n M_{n-1}+ 2 M_n + (1- \mu_n) M_n=6 f(x_{n-1},x_n,x_{n+1})$
\end{itemize}
wobei $x_{n+1} := x_n+x_1-x_0$ und $f_{n+1} := f_n$.
Mit dem Beispiel \begin{figure}

\center

\begin{tabular}{|c||c|c|c|c|}

$x_i$ & -1 & 0 & 1 & 2 \\ 
\hline 
$f_i$ & 5 & -2 & 9 & -4 \\ 

\end{tabular}\\  
\end{figure} ergibt sich daraus:
\begin{align*}
f(x_0,x_1,x_2)=9\\
f(x_1,x_2,x_3)=-12\\
\mu_i = \frac{1}{2}\\
\text{natürlicher spline:} M_0=M_3=0
\rightarrow \frac{1}{2} M_0+2M_1+\frac{1}{2} M_2 =54 \\
\frac{1}{2} M_1+ 2 M_2 + \frac{1}{2} M_3 =-72 \\
\end{align*}
Dieses Gleichungssystem ist dann zu lösen nachdem sich $C_i$ und $D_i$ durch einsetzen ergeben.
Für $n>3$: lineares Gleichungssystem für $M_n$!
Daraus resultieren 3 Splines für die 3 Intervalle:
\begin{align*}
[-1,0]:s(x)=38,4 \frac{(x+1)^3}{6} + C_i(x+0,5)+D_1\\
[0,1]:s(x)=38,4 \frac{(x-1)^3}{6} - 45,6 \frac{x^3}{6}+C_2 (x-0,5)+D_2\\
[1,2]: s(x)=-45,6 \frac{(2-x)^3}{6}+C_3 (x-1,5)+ D_3
\end{align*}
Zurück zum Beispiel $f(x)=\frac{1}{1+x^2}$. (Problem der Überschwinger).
Der kubische Spline erzeugt hierbei keine Überschwinger.

\subsection{Lineare Ausgleichsrechnung}
gegeben: $n$ Wertepaare $(x_i,y_i), i=1,...,n$ mit $x_i \neq x_j$ für $i \neq j$ \\
gesucht: stetige Funktion so dass möglichste genau $f(x_i) \approx y_i$ für alle $i$.
Wenn man die zulässigen Funktionen nicht einschränkt dann gelangt man zum Interpolationsproblem mit $f(x_i)=y_i$.
Wir schränken also den absoluten Funktionsraum ein.\\
\textbf{Def.:} gegebene Menge $F$ von stetigen Funktionen sowie n Wertepaare $(x_i, y_i)$; Funktionen $f \in F$ heißt Ausgleichsfunktion falls das Fehlerfunktional:
\begin{align*}
E(f)= \sum_{i=1}^n (y_i-f(x_i))^2
\end{align*}
minimiert wird.\\
Bemerkung: \begin{itemize}
\item häufig mehr Funktionswerte als Anpassparameter
\item Funktion f ist optimal bzgl.der 2er-Norm
\item es ist möglich die Wertepaare zu gewichten: $E(f)=  \sum_{i=1}^n (\omega_i(y_i-f(x_i))^2$ mit $\omega_i \in [0,1]$
\item für $F=$Menge aller Geraden $\rightarrow$ $f=$ Ausgleichs- oder Regressionsgerade.
\end{itemize}
\textbf{Beispiel:}
Wertepaare: \begin{figure}
\center

\begin{tabular}{|c||c|c|c|c|}

$x_i$ & 1 & 2 & 3 & 4 \\ 
\hline 
$y_i$ & 6 & 6,8 & 10 & 10,5 \\ 

\end{tabular}\\  
\end{figure}
\\
Wie lautet die Ausgleichsgerade? \\
Funktion $y=a x +b$ oder allgemein: $F:=\{a_1 f_1+a2 f_2\|a_1,a_2 \in \mathbb{R}$ mit Ansatzfunktionen $f_1(x)=x$ und $f_2(x)=1$.\\
Aufgrund der Linearität in den Koeffizienten $a_1$ und $a_2$ wird dies lineares Problem genannt.
Daraus ergibt sich das Fehlerintervall: $E(a,b)=\sum_i^4(y_i-(a x_i+b))^2$:
\begin{align*}
\frac{\partial E}{\partial a} = -2 \sum_{i=1}^n (y_i(a x_i+b))x_i=0 \\
\rightarrow a \sum_i x_i^2+ b \sum_i x_i=\sum_i x_i y_i \\
\frac{\partial E}{\partial b} = -2 \sum_{i=1}^n (y_i-(a x_i + b_i))=0 \\
\rightarrow a \sum_i x_i + b \sum_i 1=\sum_i y_i
\end{align*}
Also ein lineares Gleichungssystem für $a$ und $b$:
\begin{align*}
\begin{pmatrix}
\sum x_i^2 & \sum_i x_i \\
\sum x_i & n
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix} = \begin{pmatrix}
\sum_i x_i y_i \\
\sum_i y_i
\end{pmatrix} \\
\rightarrow
\begin{pmatrix}
30 & 10 \\
10 & 4
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix} = \begin{pmatrix}
91,6 \\
33,3
\end{pmatrix}
\end{align*}
Die allgemeine Definition des Ausgleichproblems ergibt sich wie folgt:\\
gegeben: Ansatzfunktion: $f=\sum_{i=1}^m a_i f_i$ mit $a_i \in \mathbb{R}$\\
und Wertepaare $(x_i, y_i)$ mit $i=1...n$ \\
Fehlerintervall: $E(a_n)=\sum_{i=1}^n (y_i-f(x_i))^2=\sum_{i=1}^n(y_i-\sum_{j=1}^m a_j f_j(x_i))^2=||\boldsymbol{y}-\boldsymbol{A \boldsymbol{a}}||_2$, mit:
\begin{align*}
A=\begin{pmatrix}
f_1(x_1) & f_2(x_1) & ... & f_m(x_1)\\
. & & & .\\
. & & & .\\
. & & & .\\
f_1(x_n) & ... & &f_m(x_n)
\end{pmatrix}
a = \begin{pmatrix}
a_1 \\
a_2 \\
.\\
.\\
a_m
\end{pmatrix},
y= \begin{pmatrix}
y_1 \\
.\\
.\\
.\\
a_m
\end{pmatrix}
\end{align*}

$ \boldsymbol{y} = \boldsymbol{A} \boldsymbol{a}$ nennt man das Fehlergleichungssystem und typischerweise gibt es keine Lösung, da sonst der Fehler 0 wäre $\rightarrow$ Interpolation. \\
Die Gleichungen $0=\frac{\partial E(\{a_m\})}{\partial a_i}$ mit $i=1...m$ heißen Normalgleichungen.
Das Normalgleichungssystem ergibt sich als: $A^TAa=A^Ty$\\
\textbf{Herleitung:}\\
\begin{align*}
E(\{a_m\})=(\boldsymbol{y}-\boldsymbol{A} \boldsymbol{a}))^2=(y_i-A_{ij} a_j)(y_i-A_{ik} a_k)\\
=y_i^2-A_{ij} a_j y_i-y_i A_{ik} a_k +A_{ij} a_j A_{ik} a_k \\
=y_i^2 -2 A_{ij} a_j y_i +A_{ij} a_j A_{ik} a_k \\
\frac{\partial E}{\partial a_l} = -2 A_{il} y_i + A_{il} A_{ik} a_k + A_{ij} a_j A_{il}\\
=-2 A_{il} y_i +2 A_{il} A_{ik} a_k )=0 \rightarrow A_{il} A_{ik} a_k= A_{il} y_i \\
\leftrightarrow A_{li}^T A_{ik} a_k =A_{li}^T y_i \leftrightarrow \underbrace{\boldsymbol{A}^T \boldsymbol{A}}_\text{mxm} \underbrace{\boldsymbol{a}}_\text{m Komponenten} = \underbrace{\boldsymbol{A}^T}_\text{m Komponenten}
\underbrace{\boldsymbol{y}}_\text{n Komponenten}
\end{align*}
\textbf{Bemerkungen:} Das Fehlergleichungssystem besteht aus n-Gleichungen(=Wertepaare), allerdings gibt es nur m Unbekannte(Anpassungsparamter).
$n>m \rightarrow$ das Fehlergleichungssystem ist überbestimmt, es gibt also im allgemeinen keine Lösung.\\
Man bezeichnet $r=A a -y$, als Residuumsvektor welcher im Allgemeinen nicht verschwindet.\\
Im Fall $r=0$ ist $a$ Lösung des Fehlergleichungssystems und es liegt eine perfekt Modellanpassung vor (Interpolation).
Die ist typischerweise der Fall für $m=n$ (Außer wenn Anpassungsfunktionen oder Wertepaare redundant sind...).

Also bedeutet die Lösung des Ausgleichproblems, dass das Fehlerfunktional bzgl. $\boldsymbol{a}$ minimal ist.
Die entspricht dann der Lösung des Normalengleichungssystems.
Dies entspricht einem Gleichungssystem aus $A^T A$ und $A^T y \rightarrow$ Cholesky Dekomposition möglich, der Lösungsvektor ist dann $\boldsymbol{a}$. 
Die Symmetrie für die CHolesky Dekomposition ist gegeben weil $B_{ik}=A_{ij}^TA_{jk}=A_{ki}^T A_{ji}=A_{jk}A_{ij}^T=A_{ij}^T A_{jk}$
In unserem Beispiel ergab sich: \begin{figure}[h]
\center

\begin{tabular}{|c||c|c|c|c|}

$x_i$ & 1 & 2 & 3 & 4 \\ 
\hline 
$y_i$ & 6 & 6,8 & 10 & 10,5 \\ 

\end{tabular}
\end{figure}\\
Als Anfangsfunktionen nutzen wir: $f_1(x)=x$ und $f_2(x)=1$.
\begin{align*}
A=\begin{pmatrix}
f_1(x_1) & f_2(x_1) \\
f_1(x_2) & f_2(x_2) \\
f_1(x_3) & f_2(x_3) \\
f_1(x_4) & f_2(x_4) 
\end{pmatrix}= \begin{pmatrix}
1 & 1\\
2 & 1\\
3 & 1\\
4 & 1
\end{pmatrix}, \quad A^T A = \begin{pmatrix}
1 & 2 & 3 & 4\\
1 & 1 & 1 & 1
\end{pmatrix} \begin{pmatrix}
1 & 1\\
2 & 1\\
3 & 1\\
4 & 1
\end{pmatrix}
=\begin{pmatrix}
30 & 10 \\
10 & 4
\end{pmatrix} \\
A^T y = \begin{pmatrix}
1 & 2 & 3 & 4\\
1 & 1 & 1 & 1
\end{pmatrix} \begin{pmatrix}
6 \\
6,8 \\
10\\
10,5
\end{pmatrix}= \begin{pmatrix}
91,6 \\
33,3
\end{pmatrix}
\end{align*}  
Also:
\begin{align*}
A^T A A = A^T y \leftrightarrow \begin{pmatrix}
30 & 10 \\
10 & 4
\end{pmatrix} \begin{pmatrix}
a_1 \\ 
a_2
\end{pmatrix}= \begin{pmatrix}
91,6 \\
33,3
\end{pmatrix}
\end{align*}
Also das gleiche Gleichungssystem wie vorher.
Auch bei einem nicht-linearen Funktionsansatz $f_1(x)=x^2$ und $f_2(x)=1$ ergibt sich eine lineare Ausgleichsrechnung:
\begin{align*}
A=\begin{pmatrix}
1 & 1\\
4 & 1\\
9 & 1\\
16 & 1
\end{pmatrix}, \quad A^T A = \begin{pmatrix}
1 & 4 & 9 & 16\\
1 & 1 & 1 & 1
\end{pmatrix} \begin{pmatrix}
1 & 1\\
4 & 1\\
9 & 1\\
16 & 1
\end{pmatrix}
=\begin{pmatrix}
354 & 30 \\
30 & 4
\end{pmatrix} \\
A^T y = \begin{pmatrix}
1 & 4 & 9 & 16\\
1 & 1 & 1 & 1
\end{pmatrix} \begin{pmatrix}
6 \\
6,8 \\
10\\
10,5
\end{pmatrix}= \begin{pmatrix}
291,2,6 \\
33,3
\end{pmatrix}
\end{align*}  
Wie versprochen immer noch linear:
\begin{align*}
A^T A a = A^T y \rightarrow a_1=1,296 \quad a_2=-5,59 \\
\text{mit:} \ f(x)=a_1 x^2 + a_2
\end{align*}

Für unser Beispiel ist der lineare Fit besser, also $|r|$ ist kleiner.

\subsection{Nicht-lineare Ausgleichsprobleme, Gauß-Newton Verfahren}
\textbf{Beispiel:} Wir betrachten wieder die Wertepaare $x_i$ und $y_i$ und die Funktion $f(x)=a e^{bx}$ mit den Anpassungsparametern $a$ und $b$. $f(x)$ hängt also nichtlinear von $b$ ab! \\
\textbf{Trick}: betrachte $ln(f(x))=ln(a)+b x$ und die logarithmischen Funktionswerte in der Form: \begin{figure}[h]
\center

\begin{tabular}{|c||c|c|c|c|}

$x_i$ & 1 & 2 & 3 & 4 \\ 
\hline 
$ln(y_i)$ & 1,72 & 1,52 & 2,30 & 2,35\\ 

\end{tabular}
\end{figure}\\
Damit ergibt sich: $f(x)= a x^b \rightarrow ln(f) c+ b x$

Änhlich kann verwendet werden: $f(x) = a x^b \rightarrow ln(f) = ln(a) + b ln(x) = c + b \tilde{x}$

\textbf{Echtes Beispiel:} $f(x)=a ln(x+b)$

\textbf{Definition eines allgemeinen Ausgleichssystems:} Gegeben ist eine Ansatzfunktion $f_a(a_1,... a_m, x)$ und Wertepaare $(x_i,y_i)$ mit $i=1...n$ mit dem Fehlerfunktional $E(a_I,...a_m) = \sum_{i_1}^n (y_i -f_a(a_{ij},...a_{mj},x_i))^2$\\
Dasraus ergibt sich das Ausgleichsproblem: Minimiere E im erlaubten Parameterraum.

\textbf{Bermerkung:} im allgemeinen Fall ist ein lineares Problem enthalten. Im Prinzip kann man das nicht-lineare Gleichungssystem $\frac{\partial E}{\partial a_i}=0$ mit dem Newton-Verfahren lösen; Nachteil: instabil und zweite partielle Ableitungen werden benötigt.\\
$\rightarrow$ Gauß-Newton-Verfahren \\
Definiere Vektor $\boldsymbol{g}(a_i,... a_n) = \begin{pmatrix}
y_1-f_a(a_1,... a_m, x_1)\\
.\\
.\\
.\\
y_n-f_a(a_1,... a_m, x_n)
\end{pmatrix}$ mit $f: \mathbb{R}^m \rightarrow \mathbb{R}^n$.
Das Fehlerfunktional ist dann: $E(\boldsymbol{a})=g^2, \ E:\mathbb{R}^m \rightarrow \mathbb{R}$\\
\textbf{Problem:} Finde $\boldsymbol{a}$ so dass $E(\boldsymbol{a})$ minimal wird. Dies beschreibt ein Quadramittelproblem. \\
\textbf{Idee:} $g(\boldsymbol{a})$ in $E$ wird durch eine Linearisierung ersetzt:
hierzu wird der lineare Ausdruck minimiert.
Dieser Prozess wird dann iteriert wie beim Newton-Verfahren.

Starte mit Schätzwert: $\boldsymbol{a}^{(0)}$, das lineare Ausgleichproblem ist dann definiert durch:
\begin{align*}
E(\boldsymbol{a}) = (g(\boldsymbol{a})+ D(g(a^0))(\boldsymbol{a}-\boldsymbol{a^0}))^2
\end{align*}
mit der Jacobi-Matrix $D$.
Berechne dann $\delta^{(n)}=a^{(n+1)}-a^{(n)}$ als Lösung des linearen Ausgleichsproblems, also:
\begin{align*}
\text{min}_{\delta^{(n)}}(g(a^{(n)})+D(g(a^{(n)})) \delta^{(n)})^2
\end{align*}
Setze dann $a^{(n+1)}= a^{(n)}+ \delta ^{(n)}$ und löse:
\begin{align*}
D^T(g(a^{(n)})) D(g(a^{(n)})) \delta^{(n)} = - D^T(g(a^{(n)})) g(a^{(n)})
\end{align*}
und iteriere. 

In der Praxis benutzt man ein gedämpftes Gauß-Newton-Verfahren:
\begin{itemize}
\item Berechne wieder $\delta^{(n)}$
\item Bestimmung des Dämpfungsfaktor: \\
Wähle die größte Zahl $t \in \{1, \frac{1}{2}, \frac{1}{4},...\}$ für die mit $\phi(t)=(g(a^{(n)}+t \delta^{(n)}))^2$ gilt: $\phi(t) < \phi(0)$.
\item Setze $a^{(n+1)} = a^{(n)} t \delta^{(n)}$ und iteriere.
\end{itemize}
\textbf{Bemerkung:} Der Dämpfungsfaktor $t$ erhöht den EInzugsbereich des Verfahrens, d.h. das gedämpfte Verfahren konvergiert für einen Bereich von Startvektoren $\boldsymbol{a^{(0)}}$ (vergleiche mit $t=1$, ungedämpft).

\begin{itemize}
\item Fehler wird geringer in jedem Schritt: $(\boldsymbol{g}(\boldsymbol{a}^{(n+1)}))^2=\phi(t)<\phi(0)=(\boldsymbol{g}(\boldsymbol{a}^{(n)}))^2)^2$\\
(aber: Konvergenz zum echten Minimum nicht garantiert.)
\item Abbruchbedingung: $\sqrt{(t \boldsymbol{\delta}^{(n)})^2} < \text{TOL}$\\
(aber: bedeutet nicht dass man sich im Abstand TOL vom gesuchten Minimum befinde!)
\end{itemize} 
Beispiel: \begin{figure}[h]
\center
\begin{tabular}{|c||c|c|c|c|c|}

$x_i$ & 0 & 1 & 2 & 3 & 4 \\ 
\hline 
$ln(y_i)$ & 3 & 1 & 0,5 & 0,2 &0,05\\ 

\end{tabular}
\end{figure}
mit der Ansatzfunktion: $f(x)= a_1 e^{a_2 x}$.\\
Die Verteilungsfunktion lautet also: $g(a_1,a_2)_i=y_i-a_1 e^{(a_2 x_i)}$ und die Jacobi Matrix: $D(g)= \frac{\partial g_i(a_m)}{\partial a_j}=(-e^{(a_2 x_i)}, -a_i x_i e^{(a_2 x_i)})$. Wir erwarten: $a_1>0, a_2<0$.\\
Ungedämpft $t=1$:\\
\begin{figure}[h]
\center
\begin{tabular}{|c||c|c|c|c|c|}

$i$ & 0 & 1 & 2 & 3 & 4 \\ 
\hline 
$a^{(i)})$ & 1 & 2,99 & 1,26 & 2,91 &2,98\\  
 & -1,5 & 0,32 & 0,279 & -0,86 &-1,00\\ 

\end{tabular}
\end{figure}\\
Für den Startwert $a^{(0)} = \begin{pmatrix}
2 \\
2
\end{pmatrix}$ keine Konvergenz. \\
Gedämpftes Verfahren: \\
\begin{figure}[h]
\center
\begin{tabular}{|c||c|c|c|c|}

$i$ & 0 & 1 & 2 & 3\\ 
\hline 
$a^{(i)})$ & 1 & 1,99 & 2,92 & 2,98\\  
 & -1,5 & -0,56 & -0,951 & -0,999\\ 

\end{tabular}
\end{figure}
\begin{figure}[h]
\center
\begin{tabular}{|c||c|c|c|c|c|c|}

$i$ & 0 & 1 & 2 & 3 & ... & 10 \\ 
\hline 
$a^{(i)}$ & 2 & 0,003 & 0,004 & 0,207 & ... & 10\\  
 & 2 & 2,00 & 1,75 & -0,79 & ... & -1,00\\ 

\end{tabular}
\end{figure}\\
Das gedämpfte Verfahren ist schneller und konvergenter.


\section{Numerische Differentiation/Integration}
\subsection{Numerische Differentiation}
In vielen Anwendungen $\rightarrow f'(x_0), f''(x_0)$ eienr Funktion $f(x)$ werden benötigt, aber diese Funktionen stehen nicht zur Verfügung $\rightarrow$ Näherungen!\\
\textbf{Definition:} Ableitung als Grenzwert eines Differenzenquotienten: $f'(x_0)= lim_{x \rightarrow x_0} \frac{f(x)-f(x_0) }{x-x_0} \approx \frac{f(x)-f(x_0) }{h}:= D_1 f(x_0,h)$\\
$D_1 f(x_0,h)$ ist dabei die Differenzenformel oder finite Differenz erster Ordnung.\\Herleitung genauer Differenzenformel und Fehlerschätzung mit Taylorentwicklung:
\begin{align}
f(x_0+h)=f(x_0)+ h f'(x_0)+ \frac{1}{2} h^2 f''(x_0)+...\\
f(x_0-h)=f(x_0)+ -h f'(x_0)+ \frac{1}{2} h^2 f''(x_0)+...
\end{align}
Fehlerabschätzung:\\
Aus der ersten Gleichung ist ersichtlich: $h f'(x_0)= f(x_0+h)-f(x_0)-\frac{1}{2} h^2 f''(x_0)+...$\\
\begin{align*}
\rightarrow f'(x_0)= \frac{f(x_0+h)-f(x_0)}{h}- \frac{1}{2} h f''(x_0)+... \\
f'(x_0)= D_1(f(x_0,h))- \frac{1}{2} h f''(x_0)
\end{align*}
\textbf{Definition Diskretisierungsfehler:} Bei Differenzenformeln für $f'(x_0)$ bezeichnet man den Fehler $|D_f(x_0,h)-f'(x_0)|$ als Diskretisierungs- oder Abschneidefehler. Die Formel $D_f$ hat die Fehlerordnung $k$ falls es ein $C>0$ gibt so dass für kleines $h$ gilt:
\begin{align*}
|D_f(x_0,h)-f'(x_0)| \leq C h^k
\end{align*}
Als Kurzschreibweise nutzt man für den Fehler: $\sigma(h^k)$.
In unserem Beispiel erhalten wir also: $D_1 \rightarrow |D(f(x_0,h))-f'(x_0)|\approx \frac{1}{2} h f''(x_0) \leq C h^k$ mit $k=1$.
In unserem Beispiel ist die Fehlerordnung also $k=1$! Formeln mit höherer Ordnung sind besser, da bei Halbierung von $h$ der Fehler um den Faktor $(\frac{1}{2})^k$ kleiner wird(Vorfaktoren $C$ sind also typischerweise nicht so wichtig).\\

\textbf{Herleitung einer besseren Formel für Differentiale:}\\
Berechne:
\begin{align*}
f(x_0+h)-f(x_0-h)= 2 h f'(x_0) +\frac{1}{3} h^3 f''' (x_0) + ... \\
\rightarrow D_2 f(x_0, h) := \frac{f(x_0+h)-f(x_0-h)}{2 h} \\
\text{und} D_2 f(x_0,h)-f'(x_0) = \frac{1}{6} h^2 f'''(x_0)
\end{align*}
$D_2 f$ ist also eine Differenzenformel der Fehlerordnung $k=2$ für $f'(x)$ (außerdem eine symmetrische zentrale Differenzenformel).

\textbf{höhere Ableitungen:}
Addition der Taylor-Ausdrücke:
\begin{align*}
f(x_0+h)+f(x_0-h)= 2 f(x_0)+ h^2 f''(x_0) + \frac{1}{12} h^4 f^{\text{IV}}(x_0)+... \\
\rightarrow f''(x_0)= \frac{f(x_0+h)-2 f(x_0)+f(x_0-h)}{h^2} - \frac{1}{12} h^2 f^{\text{IV}}(x_0)
\end{align*}

\textbf{Balance von Diskretisierungs- und Rundungsfehlern:}
\textbf{Beispiel:} Formel $D,f=\frac{f(x_0+h)-f(x_0)}{h}$ für $f(x)= sin(x)$ und $x_0=1$; $f'(x)=cos(x), f'(1)=0,54...$: \\

\begin{figure}[h]
\center
\begin{tabular}{c|c}

$h$ & $|D_1 f(x_0=1,h)-cos(x_0=1)|$\\ 
\hline 
$10^{-1}$ & $0,043$ \\  
$10^{-2}$ & $0,004$ \\
$10^{-3}$ & $0,0004$ \\
$10^{-4}$ & $0,00004$ \\
$10^{-5}$ & $0,000002$ \\
$10^{-6}$ & $0,000002$ \\ 
$10^{-7}$ & $0,0003$ \\
$10^{-9}$ & $0,04$ \\
$10^{-10}$ & $0,54$ \\
\end{tabular}
\end{figure}
\newpage
Daraus folgt:
$h=10^{-12} \rightarrow 0,54$\\
Für abnehmende $h$ wird der Fehler zunächst kleiner, steigt dann aber wieder an und bleibt für $h<10^{-10}= \epsilon$ konstant, da $x_0=1=1+h$ und damit $D_1(f(x_0=1, h<10^{-10}))=0$
Dilemma: Für große $h$ ergibt sich ein Diskretisierungsfehler und für kleine ein Rundungsfehler. 



\end{document}