num2.tex

\documentclass{article}

\usepackage[paper=a4paper,left=40mm,right=40mm,top=25mm,bottom=25mm]{geometry}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{tabularx}
\usepackage{paralist}
\usepackage{hyperref} %klickbares Inhaltsverzeichnis
\usepackage[rgb]{xcolor}
\usepackage{shadethm}
\usepackage{pdfpages}

\mathchardef\ordinarycolon\mathcode`\:
\mathcode`\:=\string"8000
\begingroup \catcode`\:=\active
  \gdef:{\mathrel{\mathop\ordinarycolon}}
\endgroup

\definecolor{darkgreen}{rgb}{0.0, 0.5, 0.0}

%\newcommand{\qed}{\hfill $\Box$}
\newcommand{\real}{\textrm{Re}}
\newcommand{\limn}{\lim_{n \to \infty}}
\newcommand{\einhalb}{\frac{1}{2}}
\newcommand{\teinhalb}{\tfrac{1}{2}}
\newcommand{\halbe}[1]{\frac{#1}{2}}
\newcommand{\norm}[1]{\left\Vert #1 \right \Vert}
\newcommand{\scal}[2]{\langle #1, #2 \rangle}
\newcommand{\bignorm}[1]{\big\lVert #1 \big\rVert}
\newcommand{\Bignorm}[1]{\Big\lVert #1 \Big\rVert}
\newcommand{\biggnorm}[1]{\bigg\lVert #1 \bigg\rVert}
\newcommand{\folge}[1]{({#1}_n)_{n \in \mathbb{N}}}
\newcommand{\einsdurch}[1]{\frac{1}{#1}}
\newcommand{\teinsdurch}[1]{\tfrac{1}{#1}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\ohnenull}{\setminus\{0\}}
\newcommand{\links}{\glqq$\Leftarrow$\grqq}
\newcommand{\und}{~ \text{and} ~}
\newcommand{\grad}{\text{grad}~}
\newcommand{\gdw}{\Leftrightarrow}
\newcommand{\rot}{\color{red}}
\newcommand{\gr}{\color{darkgreen}}
\newcommand{\blau}{\color{blue}}
\newcommand{\rank}{\text{rank}}
\newcommand{\ipo}{{i+1}}
\newcommand{\partiell}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\nhn}{\newline\hfill\newline}
\newcommand{\inv}{^{-1}}
\newcommand{\spann}{\textrm{span}~\!}
\newcommand{\supp}{\textrm{supp}~\!}
\newcommand{\labtag}[1]{\label{#1}\tag{#1}}
\newcommand{\ol}[1]{\overline{#1}}
\renewcommand{\l}{\left}
\renewcommand{\r}{\right}
\renewcommand{\P}{\mathbb{P}}

\renewcommand \thesection{\Roman{section}}
\renewcommand \thesubsubsection{\Roman{section}.\arabic{subsubsection}}

\makeatletter
\renewcommand*\env@matrix[1][*\c@MaxMatrixCols c]{%
  \hskip -\arraycolsep
  \let\@ifnextchar\new@ifnextchar
  \array{#1}}
\makeatother

\theoremstyle{plain}
%\newtheorem{theorem}{Theorem}[section]
%\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{definition}[theorem]{Definition}
%\theoremstyle{remark}
%\newtheorem*{remark}{Remark}
\theoremstyle{definition}
\newshadetheorem{theorem}{Theorem}[section]
\newshadetheorem{lemma}[theorem]{Lemma}
\newshadetheorem{corollary}[theorem]{Corollary}
\newshadetheorem{definition}[theorem]{Definition}
\newshadetheorem{exthm}{Excursion Theorem}
\newshadetheorem{exdef}[exthm]{Excursion Definition}
\newtheorem{motivation}[theorem]{Motivation}
\newtheorem*{remark}{Remark}
\newtheorem*{revision}{Revision}
\newtheorem*{example}{Example}

\definecolor{shadethmcolor}{HTML}{EEEEEE}

\title{Numerical Mathematics II \\ SS 2019}
\author{Lecture by Konstantin Fackeldey}

\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Basic Facts on Ordinary Differential Equations}
\begin{definition}

	An ODE of first order in some interval $I\subset \R$ is an equation of the form
\[
y'(t)=f(t,y(t)),~ t \in I
\]
where $y: I \to \C^n ,~~  y\in C^1(I)$ and $f: I\times\C^n\to\C^n$. The order is the highest derivative in the ODE. We call an ODE \textbf{explicit} if we can solve it for $y'$ and \textbf{implicit} otherwise.
\end{definition}

\begin{definition}
An ordinary differential equation of order $n$ is given as
\[
y^{(n)}(t) = f(t,y(t),y'(t),\dots ,y^{(n-1)}(t))
\]
for $t \in I\subset \R$ where $y$ is a $n$-times differentiable function on $I$ and $f:I\times(\C^n)^n\to \C^n$ is a function.

A solution $y$ of an ODE on some $J\subset I$ is a (multiple) continuously differentiable function $y: I\to \C^n$ which solves the ODE
\end{definition}
\begin{remark}
	An ODE of order $n$ can be transferred to an ODE of first order by transformation.
\end{remark}

\begin{definition}
	We call an ODE
 \[
 y^{(n)}(t)=f(t,y(t),y'(t),\dots ,y^{(n-1)}(t)), t \in I
 \]
 an \textbf{initial value problem (IVP)} for $y$ if we have additionally the constraints $y(t_0)=y_0, \dots, y^{(n-1)}(t_0)=y_{n-1}$ for $t_0 \in I$.
\end{definition}

\begin{remark}
An ODE has a swarm of solutions, IVP has specific solutions. The swarm of solutions with all constraints is called general solution.
\end{remark}


\begin{theorem}[Picard-Lindelöf]
	For $t_0\in\R,~ y_0\in\C^n,~a,b>0$ we set \[
	I=[t_0 -a, t_0+a] \und Q=\{z\in\C^n ~|~ \norm{z-y_0}_\infty \leq b\}.\] Let furthermore $F:I\times Q\to\C^n$ be continuous, with bounded components by some constant $R$ and Lipschitz-continuous in the second argument, i.e.
	\[
	\big|F_j(t,u)-F(t,v)\big| \leq L \sum_{k=1}^n |u_k - v_k|, ~~j=1,\dots,n~,~~ t\in I,~~u,v \in Q.
	\]
	Then the IVP $y'(t)=F(t,y(t)),~ y(b)=y_0$ has on $J=[t_0-\alpha, t_0 + \alpha]\subset I$ with $\alpha=\min\{a,\frac b R\}$ exactly one differentiable solution.
\end{theorem}
\begin{proof}
No Proof.
\end{proof}


\begin{remark}
The existence is local around $t_0$.
\end{remark}

\subsection{Systems of linear ODEs}

\begin{definition}
The system $y'(t)=A(t)y(t)+f(t)$ for some interval $I\subset\R$ with $A(t)=(a_{ij}(t))_{ij}\in\C^{n,n},~ a_{ij}:I\to\C \text{ for } i,j \in \{1,\dots, n\},~n\in\N, ~y: I\to \C^n$ and $f:I\to \C^n$ is a linear system of ODEs.

The function $f$ is called inhomogeneity.

The system is called \textbf{homogeneous} if $f=0$ and \textbf{inhomogeneous} otherwise.
\end{definition}

\begin{theorem}
	Let $y_1,y_n$ be two solutions of the homogeneous ODE \[
	y'(t)=A(t)y(t),~t\in I\subset\R.\] Then each linear combination \[
	\alpha y_1 + \beta y_2,~\alpha,\beta\in\C\] of $y_1$ and $y_2$ is also a solution
\end{theorem}
\begin{proof}
\[
(\alpha y_1(t)+\beta y_2(t))'=\alpha y_1'(t)+\beta y_2'(t)=\alpha A(t)y_1(t)+\beta A(t)y_2(t)=A(t)(\alpha y_1(t)+\beta y_2(t))\]
\end{proof}

\begin{remark}
The set of all solutions of homogeneous linear ODEs $y'(t)=A(t)y(t)$ form an $n$-dimensional subspace of the vector space $C^1(I,\C^n)$.

Terms of linear dependence/independence in the context of vector valued functions.
\end{remark}

\begin{definition}
Let $y_1,\dots,y_n$ be $\C^n$ valued functions. We call these functions \textbf{linear independent} on $I\subset\R$ if \[
c_1 y_1(t)+\dots+c_n y_n(t)=0 \forall t\in I\] has the only solution $c_1=\dots=c_n=0$. Otherwise we call them \textbf{linear dependent}.
\end{definition}

It is unclear, how to decide if functions are linear independent. We can make use of the so called \textit{Wronski-determinant}, which for vector valued functions $y_1,\dots y_n$ is given by
\[
W[y_1,\dots,y_n](t) := \det \l( \begin{bmatrix}
y_{11}(t) & \cdots & y_{n1}(t) \\
\vdots & & \vdots \\
y_{n1}(t) & \cdots & y_{nn}(t)
\end{bmatrix}  \r).
\]
Let us assume
\[
W[y_1,\dots,y_n](t_0)\neq 0, ~ t_0 \in I
\]
and consider the linear system
\[
\sum_{i=1}^{n }c_i y_i(t_0) = 0 \labtag{I.1}
\]
with unknown $c_1,\dots,c_n$. Due to our assumptions, $c_1=\dots=c_n=0$ is the only solution to \eqref{I.1}. Thus, $c_1=\dots,c_n=0$ is also the only solution of
\[
\sum_{i=1}^{n}c_i y_i(t) = 0 ~\forall t\in I.
\]
We have thus shown that the functions $y_1,\dots,y_n$ are also linear independent on the whole interval $I$. Hence, we have proven the following theorem.

\begin{theorem}
Let $y_1,\dots,y_n$ be linear independent solutions of $y'(t)=f(t,y(t))$. Then $W[y_1,\dots,y_n](t) \neq 0 ~\forall t\in I$.
\end{theorem}

\begin{definition}
A set of linear independent $y_1,\dots,y_n$ of $y'(t) = A(t) y(t)$ forms a \textbf{fundamental system}. The matrix
\[
Y(t) = \begin{bmatrix}
y_{11}(t) & \cdots & y_{1n}(t) \\
\vdots & & \vdots \\
y_{nn}(t) & \cdots & y_{nn}(t)
\end{bmatrix}
\]
is called \textbf{fundamental matrix}.
\end{definition}

\begin{example}
The ODE
\[
y'(t) = \begin{bmatrix} 1 & -1 \\ -2 & 0 \end{bmatrix} y(t)
\]
has the solutions
\[
y_1(t) = \begin{bmatrix} e^{2t} \\ - e^{2t} \end{bmatrix} \und y_2(t) = \begin{bmatrix} e^{-1} \\ 2 e^{-t} \end{bmatrix}.
\]
In this case, we have
\[
W[y_1,y_2](t) = \det \l( \begin{bmatrix} e^{2t} & e^{-t} \\ -e^{2t} & 2 e^{-t} \end{bmatrix}  \r) = 3 e^t \neq 0.
\]
\end{example}

\begin{theorem}
Each linear differential system $y'(t) = A(t) y(t)$ has a fundamental system. The general solution of the system is given as
\[
y= c_1 y_1 + \dots + c_n y_n = Yc,
\]
where $Y=[y_1,\dots,y_n]$ is the fundamental matrix and $c\in \C^n$.
\end{theorem}
\begin{proof}
No proof.
\end{proof}

\begin{theorem}[Picard-Lindel\"{o}f for linear systems]
Let $A$ and $f$ be continuous on $I\subseteq \R$. Then, the IVP $y'(t) = A(t) y(t) + f(t),~ y(t_0) = y_0$ has a unique solution.
\end{theorem}
\begin{proof}
We only give a sketch of the proof. Since for a solution $y$ it holds
\[
y'(t) = A(t) y(t) + f(t),
\]
by integration we obtain
\[
\int_{t_0}^{t}y'(s)~ ds = \int_{t_0 }^{t }A(s) y(s) + f(s) ~ds.
\]
Hence,
\[
y(t) = y(t_0) + \int_{t_0 }^{t }A(s) y(s) + f(s) ~ ds.
\]
Now, by a fixpoint iteration, we get the \textit{Picard iteration} by
\begin{align*}
y_0(t) & = y(t_0) \\
y_1(t) & = y(t_0) + \int_{t_0 }^{t }A(s) y_0(s) + f(s) ~ds \\
& \vdots \\
y_{n+1}(t) & = y_(t_0) + \int_{t_0 }^{t }A(s) y_n(s) + f(s) ~ ds.
\end{align*}
It can be shown that this sequence converges to the solution pointwise.
\end{proof}

\begin{remark}
For linear ODEs, the existence and uniqueness of the solution is not local but global. Consider the special solution $y_p$ of
\[
y'(t) = A(t) y(t) + f(t) \labtag{I.2}
\]
and the homogeneous solution $y_h$ of
\[
y'(t) = A(t) y(t). \labtag{I.3}
\]
Then we can write all the solutions $y$ of \eqref{I.2} as
\[
y = y_h + y_p = Yc + y_p.
\]
The special solution $y_p$ can be found using \textit{variation of constants}, i.e. choosing $c=c(t)$. {\rot The ansatz is missing.}
\end{remark}

\begin{theorem}
The IVP
\begin{align*}
y'(t) & = A(t) y(t) + f(t) \\
y(t_0) & = y_0
\end{align*}
has the unique solution
\[
y(t) = Y(t) \cdot Y\inv(t_0) y_0 + Y(t) \cdot \int_{t_0 }^{t } Y^{-1}(s)f(s)~ds,
\]
where $Y(t)$ is the fundamental solution of the corresponding homogeneous system.
\end{theorem}
\begin{proof}
No proof.
\end{proof}

\subsection{Separated ODEs}
\begin{definition}
A ODE of the form
\[
y'(t) = g(y(t)) h(t),~ t \in I\subseteq \R,
\]
where $g,h$ are continuous functions is called \textbf{separable}.
\end{definition}

\begin{remark}
The function $g(y)$ depends only on $y$, whereas $h(t)$ depends only on $t$. We can thus integrate separatly. This gives the so called \textit{separation of variables}. In an autonomous case, i.e. $f(y,t) = f(y)$, we have no dependence on $t$ and hence $h(t) = 1$ for all $t$. In this sense, an autonomous system is always separable.
\end{remark}

\setcounter{subsection}{2}
\subsection{Qualitative Behaviour of ODEs}
\begin{example}
Let us consider the $n$-dimensional non autonomous system of first order
\begin{align*}
	y'(t)&=f(t,y(t))\\
	y(t_0)&=y_0
\end{align*}
where $f: D \to \R^n, ~D\subset I\times \R^n,~ t_0 \in I,~ I\subset \R$. The questions we are dealing with are:
\begin{enumerate}
\item Why only first order?
\item What is the relation between a non-autonomous and an autonomous system?
\end{enumerate}
\end{example}

The reason behind 1. is that any ODE of $n$-th order can be transformed into a $n$-dimensional ODE of first order. Consider the ODE
\[
x^{(n)}=F(t,x(t),x'(t),\dots,x^{(n-1)}(t))
\]
and define a vector $y$ with its components $y_i, ~ i=1,\dots,n$ by
\[
y_i(t)=x^{(i-1)}(t)
\]
and a vector field $f(t,y)$ by
\[
f(t,y) = \Big(t,y_1,y_2,\dots,y_n,F(t,y_1,y_2,\dots,y_n)\Big)^T
\]
Then the ODE of $n$-th order is equivalent to $y'(t)=f(t,y)$.\newline

A system of the form $y'(t)=f(t,y)$ is called an \textbf{non-autonomous system}, a system of the form $y=f(y)$ is called \textbf{autonomous}. We can transform a non-autonomous system to an autonomous system.

Consider the ODE
\[
y'(t)=f(t,y) \und y(t_0) = y_0.
\]
We set
\[
z=\begin{bmatrix} y \\ s \end{bmatrix} \und \hat{f}=\begin{bmatrix} f(s,y) \\ 1 \end{bmatrix}, ~~ s\in \R
\]
Then
\[
z'(t)=\hat{f}(z(t)), ~ z(t_0)=z_0=\begin{bmatrix} y_0\\t_0 \end{bmatrix}
\]
is an autonomous system.

In short, each ODE in $\R^n$ can be transformed to an autonomous ODE in $\R^{n+1}$

\begin{remark}
In the theorem of Picard-Lindelöf the ODE is of the form $f(t,y)$, where $y$ has to be Lipschitz-continuous.

In the autonomous system the right hand side looks like $f(y(t))$, where $t$ and $y$ have to be Lipschitz-continuous.
\end{remark}

\subsubsection*{Analytic Continuation}
"Local solutions can be spread onto a maximum time interval."\newline

\setcounter{theorem}{13}
\begin{definition}[Local Lipschitz]
A function $f: X \to Y$ is local Lipschitz in $x\in X$ if the exists a neighbourhood $U_x \subseteq  X$ around $x$ such that $f|_{U_x}$ is Lipschitz-continuous.
\end{definition}

For $G:=I\times Q$ with $I=[t_0-a, t_0+a], ~ Q=\{z\in\C~|~ \norm{z-y_0}\leq b\}$ with $a,b>0$ the theorem of Picard-Lindelöf gives for local Lipschitz $f$ the existence of a solution $y_0(t)$ of the IVP
\begin{align}
	y'(t)&=f(t,y(t)) \tag{1.4}\\
	y(t_0)&=y_0\notag
\end{align}
on some (small) inverval $I_0=[t_0-a_0,t_0+a_0]$ with $a_0=a>0$.

We will have a look at what happens if we apply the theorem of Picard-Lindelöf on one side of the interval $I_0$. Let now be $t_1:=t_0+a_0$ and $y_1 = y_0(t_1)$. We then have that $(y_1,t_1)\in G$ and according to Picard-Lindelöf we know that the IVP with $y(t_1)=y_1$ has a unique solution $y_1(t)$ on $I_1:=[t_1-a_0,t_1+a_1]$ where $a_1>0$.

Due to the uniqueness of the solution if hold $y_0(t)=y_1(t)$ on $I_0\cap I_1$ we are defining a continuoation of our solution on the greater interval.

It holds
\[
y_{+}(t)=y_0(t) ~\text{for}~ t\in[t_0,t_1]
\]
and
\[
y_{+}(t)=y_1(t) ~\text{for}~ t \in (t_1,t_1+a_1]
\]
analogue for $y_{-}(t)$. Thus there exists a unique solution on the interval $[t_0, t_0+a_0+a_1+\dots]$ if $\sum_{k=0}^{\infty}a_k<\infty$. If $\sum_{k=0}^{\infty}a_k$ diverges, the solution exists globally in forward time.

\begin{remark}
It can happen that $a_n$ can arbitrary small when $(t_k,y_{+}(t_k))$ approaches the boundary of $G$. Then either $\norm{f((t_k),y_{+}(t_k))}$ or the Lipschitz-constant $L$ might get arbitrary large.
\end{remark}

\begin{definition}
Let $f:G \to \R^n$ be continuous and local Lipschitz with respect to $y$ and let $(t_0,y_0)\in G$. Let furthermore $t_{\pm}:=t_{\pm}(t_0,y_0)\in\R$ be defined as
\begin{align*}
t_{+}&=\sup \{\tau>t_0 ~|~ \text{there exists a continuation}~ y_+ \text{ of (1.4) on } [t_0,\tau]\}\\
t_{-}&=~\inf \{\tau>t_0 ~|~ \text{there exists a continuation}~ y_- \text{ of (1.4) on } [t_0,\tau]\}.
\end{align*}
The interval $(t_-,t_+)$ is the largest interval of existence of the IVP with some initial point $y(t_0)=y_0$.

The \textbf{maximum solution} $y(t)$ is
\[
y(t)=\begin{cases} y_+(t) \text{ for } t\in[t_0,t_+) \\ y_-(t) \text{ for } t\in (t_-,t_0].\end{cases}
\]
\end{definition}

\begin{example}
Consider
\[
y'=y^2, ~ y(t_0)=y(0)=1, ~ y(t)=\frac{1}{1-t}.
\]
Then we have $(t_-,t_+)=(-\infty, 1) \text{ or } (1,\infty)$.
\end{example}

\begin{remark}
In case of $t_+<\infty$ the maximum solution approches for $t\to t_+$, it can then happen that $\norm{y(t)}$ is unbounded. This is also called "blow up".
\end{remark}

\subsubsection*{Solutions and Initial Data}
"What is the influence of a perturbation in $f,~y_0$ or $t_0$ on the solution?"\newline

To consider this, we need the following Lemma.
\begin{lemma}[Grönwall-Lemma]
Let $I=[a,b]\subseteq \R$ and $g:I\to\R$ be a continuous function. If
\[
0\leq g(t)\leq \delta+\gamma \int_{a}^{t }g(x)~dx
\]
holds for all $t\in I, ~\delta,\gamma >0$, then it holds
\[
g(t)\leq\delta e^{\gamma(t-a)}.
\]
\end{lemma}
\begin{proof}
We set
\[
\varphi(t)=\delta+\gamma \int_{a}^{t }g(x)~dx.
\]
Then we have
\[
\varphi'(t)=\gamma\cdot g(t) \leq \gamma\varphi(t).
\]
Since
\[
\Big(\varphi\cdot e^{-\gamma t}\Big)' = \varphi'\cdot e^{-\gamma t} + \varphi \cdot (-\gamma) e^{-\gamma t}= e^{-\gamma t}\Big(\varphi'(t)-\gamma\varphi(t)\Big)\leq 0
\]
we have that $\varphi e^{-\gamma t}$ is monotone falling. It thus follows
\[
g(t)\cdot e^{-\gamma t}\leq \varphi(t)\cdot e^{\gamma t} \leq \varphi(a)\cdot e^{-\gamma a}= \delta \cdot e^{-\gamma a}
\]
for all $t\geq a$.
\end{proof}
The Grönwall-Lemma allows us to prove the following theorem.
\begin{theorem}[Dependence on initial data]
Let $D\subset I\times \R^n$ be open, $f:D\to\R^n$ continuous and local Lipschitz with respect to $y$ and $(t_0,y_0) \in D$. If the solution of
\begin{align*}
y'(t)&=f(t,y(t))\\
y(t_0)&=y_0, ~y_0\in\R^n
\end{align*}
exists for all $t\in I=[a,b]$ then for each $\varepsilon>0$ there exists a $\delta>0$ such that \begin{enumerate}[(i)]
\item If $\norm{y_0-z_0}<\delta$ there also exists a solution of
\begin{align*}
z'(t)&=f(t,z(t))\\
z(t_0)&=z_0, ~z_0\in\R^n
\end{align*}
for $t \in I$.
\item It holds
\[
\max_{t\in I}\norm{y(t)-z(t)}<\varepsilon.
\]
\end{enumerate}
\end{theorem}
\begin{proof}
Since $D$ is open, there exists a $\bar\delta>0$ and a compact set
\[
K:=\{(t,z(t))~|~t\in I, ~ \norm{y(t)-z(t)}\leq\bar\delta\}\subset D.
\]
On $K$ the function $f$ is Lipschitz (with respect to $y$) with a Lipschitz-constant $L$. Let now $\delta<\bar\delta$ and $\norm{y_0-z_0}<\delta$. Then for all $t_0,t\in[a,b]$ it holds
\[
\norm{z(t)-y(t)}\leq \delta+L \int_{t_0}^{t }\norm{y(x)-z(x)}~dx.
\]
This can be seen by the integral representation of $y(t)$. Applying Grönwall's Lemma with $\gamma=L$ yields
\begin{equation}
\norm{y(t)-z(t)}\leq\delta\cdot e^{L(t-t_0)} \labtag{I.5}
\end{equation}
and by choosing $\delta\leq\bar\delta\cdot e^{L(a-b)}$ it holds $\norm{y(t)-z(t)}\leq \bar\delta$ for all $t\in I$. Thus it holds $(t,z(t))\in K$ for $t\in[a,b]$ and hence we have shown \textit{(i)}.

By choosing $\delta<\varepsilon\cdot e^{L(a-b)}$ it follows \textit{(ii)}.
\end{proof}
\begin{remark}
We have thus shown, that the solution $y(t)$ of the IVP with initial value $y(t_0)=y_0$ depends continuously on the initial data. The solution is often written as $y(t;t_0,y_0,f)$.
\end{remark}
\begin{example}
Let us consider the ODE
\begin{align*}
y'&=\lambda y, ~ \lambda \in \R\\
y(0)&=y_0
\end{align*}
Here we have $L=|\lambda|$. The equation \eqref{I.5} gives
\[
|y(t)-z(t)|\leq e^{|\lambda|\cdot t} |y_0-z_0|.
\]
For $\lambda<0$ we know that $|y(t)-z(t)|$ decreases exponentially.
\end{example}
\subsection{Stability and Flow}
\subsubsection*{Vector field}
A solution of an ODE is a function $y:I\to\R^n$ which is differentiable on $I$. Its graph $\{(t,y(t))~|~t\in I\}$ is a differentiable curve in $\R^{n+1}$ also known as \textit{solution curve} or \textit{integral curve}. In each point $(t,y(t))$ the direction of the tangent is given by the $(1,f(t,y(t)))$. In other words, $f$ is assigning a direction to each point.

\subsubsection*{Stability and small perturbations}
Consider
\[
y'(t)=f(t,y(t)), ~ y(t_0)=y_0.
\]
We are now interested in a comparison of different solutions for $t\in[t_0,\infty)$ with respect to the initial condition. We denote the solution by $y(t)=y(t,t_0)$.

Stability means that $y(t_0)=\tilde y$ with $\tilde y$ near by $y_0$. The question we are dealing with is "How does $y(t,\tilde y)$ behave in comparison with $y(t,y_0)$?".\newline

Let us consider an autonomous ODE, i.e. an ODE of the form $y'(t)=f(y(t))$.

\begin{definition}[Equilibrium Point]
A point $\bar y\in D\subset\R^n$ is called an equilibrium point of a mapping $f:D\to\R^n$ if $f(\bar y)=0$. The constant solution $y(t)=\bar y$ is the only solution with $y(t_0)=\bar y$.
\end{definition}
\begin{remark}
Other names for equilibrium points are fixed points, equilibria and stationary points.
\end{remark}
\begin{definition}[Stability and asymptotic stability]
An equilibrium point is \textbf{stable} (in the sense of Ljapunov) if for each $\varepsilon>0$ there exists a $\delta>0$ such that for $t\geq t_0$ and for all trajectories $y(t)$ with $\norm{y(t_0)-\bar y}\leq \delta$ it holds that
\[
\norm{y(t)-\bar y}\leq \varepsilon.
\]

An equilibrium point is \textbf{instable} if it is not stable.

An equilibrium point $\bar y$ is \textbf{asymptotic stable} if there exists a neighbourhood $U_{\bar y}$ of $\bar y$ such that
\[
y(t_0)\in U_{\bar y} \Rightarrow \lim_{t\to\infty} y(t)=\bar y.
\]
In this case $\bar y$ is called a \textbf{sink}.

An equilibrium point $\bar y$ is a \textbf{spring} if for each solution $y(t)$ with $y(t_0)\in U_{\bar y}$ and $y(t_0)\neq \bar y$ there exists a $t_1>t_0$ such that $y(t)\not \in U_{\bar y}$ for all $t\geq t_1$.
\end{definition}
\begin{example}
Consider an ODE in $\R^1$ given by $y'(t)=f(t,y(t))$. The equilibrium point is asymptotic stable if in $U_{\bar y}$ it holds that
\[
f(y)<0 \text{ for } y<\bar y ~~\und~~ f(y)>0 \text{ for } y>\bar y.
\]
\end{example}
\begin{definition}[Stability of solutions]
Let $y(t;y_0)$ be a solution of $y'(t)=f(y(t)), ~y(t_0)=y_0 ~\forall t\geq t_0$. Then the solution is \textbf{stable} if for each $\varepsilon>0$ there exists a $\delta>0$ such that
\[
\norm{y_0-\tilde y_0}\leq \delta \Rightarrow \norm{y(t;y_0)-y(t,\tilde y_0)}<\varepsilon
\]
for all $t>t_0$.
The solution is \textbf{attractive} if there exists a $\delta >0$ such that
\[
\norm{y_0-\tilde y_0}<\delta \Rightarrow \lim_{t\to\infty}\norm{y(t;y_0)-y(t,\tilde y_0)}=0.
\]
The solution is \textbf{asymptotic stable} if its stable and attractive.
\end{definition}
\subsubsection*{Flow and Dynamical System}
A Dynamical System is a mathematical model to understand a time independent (autonomous) process. This process shall not depend on the initial time but only on the initial state. Formally, a dynamical system is triple $(T,S,\Phi)$ where $T$ is the time space, $S$ is the state space and $\Phi:T\times S\to S$ is the flow. The time space can either be discrete ($T=\N$) or continuous ($T=\R,~ S=\R^n$). This dynamical system is described by an ODE: The entity of all solutions of an ODE is a dynamical system
\[
y'(t)=f(y)
\]
where $f$ is a differentiable vector field.

\begin{definition}[Flow of an autonomous ODE]
The flow $\Phi(t,y_0)$ or $\Phi_{t}(y_0)$ of an autonomous ODE
\[
y'(t)=f(y(t)), ~ y(t_0)=y_0
\]
is a mapping $\Phi:\R^{n+1}\to \R^n, ~ \Phi(t,y_0)=y(t)$ and with the following properties:
\begin{enumerate}[(i)]
\item $\Phi(t_0,y_0)=y_0$ for all $y_0\in\R^n$
\item $\Phi(t_1+t_2,\cdot) = \Phi\Big(t_2,\Phi(t_1,\cdot)\Big)$ for $t_1,t_2\in \R$.
\end{enumerate}
\end{definition}

\begin{remark}\mbox{}
\begin{itemize}
\item $\Phi(t,y_0)$ is the solution of the ODE $y'(t)=f(y(t))$ which starts in $y_0$ at $t_0$.
\item $\Phi:\R^{n+1}\to\R^n$ is differentiable, i.e. $\Phi(t,y_0)$ is a $C^1$/function and it holds
\[
\partiell{}{t}\Phi(t,y_0)=f\Big(\Phi(t,y_0)\Big).
\]
\end{itemize}
\end{remark}

\begin{example}
For the ODE
\begin{align*}
y'(t)&=Ay(t)\\
y(t_0)&=y_0
\end{align*}
with $A\in \R^{n,n}$ it holds
\[
\Phi(t,y_0)=e^{At}y_0
\]
for all $t\in \R$.
\end{example}

\begin{lemma}
Under the assumptions of the theorem of Picard-Lindelöf on the ODE
\[
y'(t)=f(y(t))
\]
the solutions $y_1$ and $y_2$ of different initial conditions do not intersect.
\end{lemma}
\begin{proof}
Let us assume towards a contradiction that we have two solutions $\Phi(t_1,y_1)$ and $\Phi(t_2,y_2)$ with different initial conditions which intersect at $y^*$, i.e.
\[
\Phi(t_1,y_1)=\Phi(t_2,y_2)=y^*.
\]
We define
\[
v(t):=\Phi(t+t_1,y_1)=\Phi\big(t,\Phi(t_1,y_1)\big)=\Phi(t,y^*)
\]
and
\[
w(t):=\Phi(t+t_2,y_2)=\Phi\big(t,\Phi(t_2,y_2)\big)=\Phi(t,y^*).
\]
Then by the theorem of Picard-Lindelöf it follows that
\[
v(t)=w(t)
\]
what ends the proof.
\end{proof}

\noindent By
\[
\mathcal{O}(y_0):=\{y\in\R^n~|~\exists t \in \R:~ y=\Phi(t,y_0)\}
\]
we denote the image of the mapping $t\to \Phi(t,y_0)$. The set $\mathcal{O}(y_0)$ is called \textbf{trajectory} or \textbf{orbit}.

\begin{example}[Predator-Prey-Model, Räuber-Beute-Modell]
Let $x$ represent the number of prey (maybe a goat) and $y$ the number of the predators (maybe a wolf). We can model
\begin{align*}
 x'&=x(a-by) \labtag{I.6} \\
 y'&=y(-c+dx)
\end{align*}
where $a,b,c,d\in\R_{>0}$. In the absence of preadotrsthe number of prey is growing exponentially. An increase in the number of predators means a decrease in the number of preys. Note that the decrease of the preys is propotional to $x\cdot y$. In the absence of preys, the predators die. An increase in the number of preys means an increase in the number of predators.

Also note that we assume that the wolf only eats goats and that no further enemies of the goat exist.

These equations belong to the Lotha-Volterra equations.

The origin $(0,0)$ is the only equilibrium point on the boundary of the state space $\R^{2}_{\geq 0}$. In the interior of $\R^2_{\geq 0}$ there exists also only one equilibrium point which is given by $(\bar x, \bar y)=(\tfrac{c}{d},\tfrac{a}{b})$.

The curves of the solutions are closed. To see this, reconsider \eqref{I.6}. Using simple calculations we get
\[
x'\Big(\frac{c }{x}-d\Big)=(a-by)(c-dx)
\]
and
\[
y'\Big(\frac{a }{y}-b\Big)=(-c+dx)(a-by).
\]
By adding up, we obtain
\[
\Big(\frac{c }{x}-d\Big)x' + \Big(\frac{a }{y} - b\Big)y' = 0
\]
or (using the method of \textit{scharf hinsehen})
\[
\partiell{}{t}\Big(c\ln(x)-dx+a\ln(y)-by\Big)=0.
\]
Setting
\[
B(x):=\bar x \cdot \ln(x) - x ~~\und~~ R(y):=\bar y \cdot \ln(y) - y
\]
it holds for $V(x,y):=dB(x)+bR(y)$ that
\[
\partiell{}{t} V\big(x(t),y(t)\big)=0
\]
or $V(x,y)$ is constant along the trajectories of the solutions. We see that $V(x,y)$ is a conserved quantity (\textit{Erhaltungsgröße}) taking its maximum in the equilibrium point $(\bar x, \bar y)$. This point is stable, too (Homework).
\end{example}

Let us now consider $V:D\to \R, ~D\subseteq \R^n$ such that in $D$ there exists a equilibrium point $\bar y$ of the system $y'=f(y)$. Taking the derivative of $V$ along the solution $y(t)$ we obtain
\[
V'(y(t))=\partiell{}{t}V(y(t))=\nabla \Big(V\cdot y'(t)\Big) = \nabla V\big(f(y(t)\big).
\]
If $V'\leq 0$, then $V$ is a monotone falling function along all solutions $y(t)\in D$.

\setcounter{theorem}{23}
\begin{theorem}[Ljapunov-Stability]
Let $\bar y\in D\subseteq \R^n$ be an equilibrium point of $y'=f(y)$. Let further $V:D\to\R$ be a differentiable function on an open set $D$ and let $V(\bar y)=0$ and $V(y)>0$ for $y\neq \bar y$ and
\[
V'=\partiell{}{t}V \leq 0 ~~\text{on}~~ D\setminus \{\bar y\}.
\]
Then the equilibrium point $\bar y$ is stable. If we have $V'<0$ then $\bar y$ is asymptotic stable.
\end{theorem}
\begin{proof}
No proof.
\end{proof}
\begin{remark}
The function $V$ from theorem I.24 is called Ljapunov-function.
\end{remark}

\newpage
\section{Numerics of ODEs}
\begin{motivation}
     In the following we only consider first order ODEs for a bounded interval $[a,b]\subseteq R$ and a given function $f:[a,b]\times \R \to \R$. We seek for a function $y: [a,b]\to \R$ such that
     \footnote {We assume in \eqref{II.1} that $f$ is sufficiently small, such that all necessery (Taylor-)expansions can be built and we also have uniqueness and existence of a solution for the IVP.}
     \begin{align*}
     y'(t)&= f(t,y(t)) ~~ \forall t\in[a,b] \labtag{II.1}
     \end{align*}
     with initial condition
     \begin{align*}
     y'(a) = \hat y. \labtag{II.2}
     \end{align*}
     We devide the interval $[a,b]$ by
     \[
     a=t_0 < t_1 < \dots < t_n=b, ~~ \Delta t_i= t_{i+1}-t_i.
     \]
     At the beginning we only consider an equidistant mesh, i.e. $\Delta t_i$ is constant. Later we also consider variable meshsizes, since there might exist solutions where variable meshsizes can be helpful.
     We write
     \[
     \Delta t = \frac{b-a}{n} \und t_i = t_0+i\cdot\Delta t.
     \]
     Given a starting value $y_0$ we compute our approximations $y_i$ of the exact solution $y(t_i)$ evaluated at $t_i$.
\end{motivation}
\subsection{Two different schemes}
\subsubsection*{Difference method}
Replace the tangent of $y$ at $t_i$ by a secant with respect to $t_i$ and $t_{i+1}$, i.e.
\[
y'(t_i)=\frac{y(t_{i+1})-y(t_i)}{\Delta t}.
\]
Inserting this into the ODE gives
\[
\frac{y(t_{i+1})-y(t_i)}{\Delta t}\approx f(t,y(t)).
\]
This leads to the \emph{explicit Euler-Method}
\[
y_{i+1}= y_i + \Delta t\cdot f(t_i,y_i), ~~ i=0,\dots,n-1.
\]
\subsubsection*{Integration method}
We are using the equation
\[
y(t_{i+1})-y(t_i)=\int_{t_i }^{t_{i+1}}y'(\tau)~d\tau=\int_{t_{i}}^{t_{i+1}}f(\tau,y(\tau))~d\tau.
\]
Applying the quadrature rule leads to
\[
\int_{t_i}^{t_{i+1}}f(\tau,y(\tau))~d\tau \approx (t_{i+1}-t_i)\cdot f(t_{i+1},y(t_{i+1})).
\]
The \emph{implicit Euler-Method} follows by that as
\[
y_{i+1}=y_i+\Delta t \cdot f(t_\ipo,y(t_{i+1})), ~~i=0,\dots,n-1.
\]
\subsection{One-Step Methods}
"For computing $y_{i+1}$ of $y$ we only use the information at $t_i$."
\begin{definition}[One-Step Method]
A method for approximating the IVP \eqref{II.1} and \eqref{II.2} of the form
\[
y_{i+1}=y_i + \Delta t ~\Phi(t_1,y_i,y_\ipo,\Delta t)
\]
with some given starting value $y_0$ at $t_0$ and an incremental function (\textit{Verfahrensfunktion})
\[
\Phi:[a,b]\times \R\times\R\times \R_+ \to \R
\]
is called a \textbf{one-step method}. We call it \textbf{explicit} if $\Phi$ depends not on $y_\ipo$ and \textbf{implicit} otherwise.
\end{definition}
\begin{example}
For the explicit Euler-Method the incremental function $\Phi$ is
\[
\Phi(t_i,y_i,y_\ipo,\Delta t)= f(t_1,y_i).
\]
For the implicit Euler-Method the incremental function $\Phi$ is
\[
\Phi(t_i,y_i,y_\ipo,\Delta t)= f(t_{i+1},y_{i+1}).
\]
\end{example}
Note that in the following we use an abuse of notation: In the explicit case we write $\Phi(t_1,y_1,\Delta t)$. \newline\hfill

\noindent But how do we measure the quality of our approximation?

\begin{definition}[local discretization error (consistency)]
A one-step method is \textbf{consistent of order $p\in\N$} if for an ODE \eqref{II.1} with some solution $y$ and the local discretization error
\[
\eta(t,\Delta t)=y(t)+\Delta t \cdot \Phi\Big(t,y(t),y(t+\Delta t),\Delta t\Big) - y(t+\Delta t)
\]
for $t\in [a,b]$ and $0\leq \Delta t \leq b-t$ it holds
\[
\eta(t,\Delta t)=O(\Delta t^{p+1}) ~\text{ as }~ \Delta t \to 0.
\]
In case of $p=1$ we say that the method is \textbf{consistent}.
\end{definition}
\begin{revision}
The Landau-Notation for functions $f$ and $g$ is defined as follows:\newline
It holds "$f(x)=O(g(x))$ for $x\to a$" if $\Big|\frac{f(x)}{g(x)}\Big|$ is bounded when $x\to a$. Furthermore it holds "$f(x)=o(g(x))$ for $x\to a$" if $\lim_{x\to a} \frac{f(x)}{g(x)}=0$. We make use of an abuse of notation by writing the equality sign, since formally $O(g(x))$ and $o(g(x))$ are sets.
\end{revision}
\begin{remark}
For a consistent method it holds
\begin{align*}
\lim_{\Delta t\to 0} \Phi\Big(t,y(t),y(t+\Delta t), \Delta t\Big) &= \underbrace{\lim_{\Delta t \to 0} \frac{\eta(t,\Delta t)}{\Delta t}}_{=0}+\lim_{\Delta t \to 0} \frac{y(t+ \Delta t)- y(t)}{\Delta t}\\&= y'(t)= f(t,y(t)).
\end{align*}
\end{remark}

\setcounter{theorem}{2}
\begin{theorem}[Consistence of the explicit Euler-Method]
The explicit Euler-Method is consistent of order $p=1$.
\end{theorem}
\begin{proof}
Expansion of $y$ in $t$ gives
\begin{align*}
y(t+\Delta t)
&=y(t)+y(t)'\cdot\Delta t + \frac{y''(\varrho)}{2}\Delta t^2,~~\varrho \in [t,t+\Delta t]\\
&=y(t)+f(t,y(t))\cdot \Delta t + \frac{y''(\varrho)}{2}\Delta t^2.
\end{align*}
It thus follows
\begin{align*}
\eta(t,\Delta t)
&=y(t)-\Delta t \cdot f(t,y(t)) - y(t+\Delta t)\\
&= -\frac{\Delta t^2}{2}y''(\varrho) = O(\Delta t^2)
\end{align*}
for $\Delta t \to 0$, since $y''$ is bounded in $[t,t+\Delta t]$.
\end{proof}

\begin{definition}[Convergence of one-step methods]
A one-step method with starting value $y_0=y(0)+O(\Delta t^p),~\Delta t \to 0$ is \textbf{convergent of order $p\in\N$} with respect to the IVP \eqref{II.1} and \eqref{II.2} if for the approximation $y_i$ of the solution $y(t_i)$ the \textbf{global approximation error}
\[
e(t_i,\Delta t)=y(t_i)-y_i
\]
for all $t_i,~i=1,\dots, n$ meets
\[
e(t_i,\Delta t)=O(\Delta t^p), ~~ \Delta t \to 0
\]
In case of $e(t,\Delta t)=O(1)$ we call the method \textbf{consistent}.
\end{definition}

\begin{remark}
Note, that in $e(t_i,\Delta t)$ all $\eta(t,\Delta t)$ are summed up.
\end{remark}

\begin{lemma}[technical Lemma]
Let $\eta_i,\varrho_i,z_i\in\R_{\geq 0}$ for $i=0,\dots,m-1$ and $z_m\in\R$ and it holds
\[
z_\ipo\leq (1+\varrho_i)z_i + \eta_i
\]
for $i=0,\dots,m-1$. Then it holds
\[
z_{i+1} \leq \l(z_0 + \sum_{k=0}^i\eta_k\r) e^{\sum_{k=0}^i\varrho_k}
\]
for $i=0,\dots,m-1$.
\end{lemma}
\begin{proof}
We prove the statement by induction on $i$. For $i=0$ the claim is true. Hence let the statement be valid for a $i-1$. Then we have
\begin{align*}
z_\ipo
& \leq (1+\varrho_i) z_i + \eta_i \\
& \leq (\underbrace{1+\varrho_i}_{\leq e^{\delta_i}}) \cdot \l(z_0 + \sum_{k=0}^{i-1}\eta_k\r)e^{\sum_{k=0}^{i-1}\varrho_k} + \eta_i \\
& \leq \l(z_0 + \sum_{k=0}^{i-1}\eta_k\r) e^{\sum_{k=0}^{i-1}\varrho_k} + \eta_i \\
& \leq \l(z_0 + \sum_{k=0}^{i} \eta_k \r)\cdot e^{\sum_{k=0}^{i}\varrho_k},
\end{align*}
what ends the proof.
\end{proof}

\begin{theorem}[Convergence of one-step methods]
Let $\Phi$ be an incremental function of a one-step method for the IVP \eqref{II.1} and \eqref{II.2} with
\begin{align*}
\l| \Phi(t,u,w,\Delta t) - \Phi(t,v,w,\Delta t)\r| & \leq L |u-v| \labtag{II.3} \\
\l| \Phi(t,w,u,\Delta t) - \Phi(t,w,v,\Delta t)\r| & \leq L |u-v| \labtag{II.4}
\end{align*}
with $L\in \R$. Then it holds for $\Delta t < \tfrac{1}{L}$ that
\[
\l| e(t_{i+1},\Delta t)\r| \leq \l( |e(t_0,\Delta t)| + \frac{(t_{i+1}-t_0)}{1-\Delta t\cdot L}\cdot \frac{\eta(\Delta t)}{\Delta t}\r) e^{2 \cdot \frac{t_{i+1}-t_0}{1-\Delta t}\cdot L}\labtag{II.5}
\]
for $i=0,\dots,n-1$, where
\[
\eta(\Delta t) := \max_{j=0,\dots,n-1} |\eta(t_j,\Delta t)|.
\]
\end{theorem}
\begin{proof}
Reconsider that
\[
\eta(t_i,\Delta t)= y(t_i)+\Delta t\Phi\l(t_i,y(t_i),y(t_i+\Delta t),\Delta t\r) - y(t_{i+1}).
\]
Rearranging gives
\[
y(t_\ipo)= y(t_i) + \Delta t\Phi\l(t_i,y(t_i),y(t_i+\Delta t),\Delta t\r) - \eta(t_i,\Delta t).
\]
Consider now
\begin{align*}
e(t_\ipo,\Delta t)
& = y(t_\ipo)-y_\ipo  \\
& = y(t_i)+\Delta t \Phi\l(t_i, y(t_i),y(t_\ipo),\Delta t\r) - \eta(t_i,\Delta t) \\
& ~~~ - y_i -\Delta t \Phi \l(t_i, y_i, y_\ipo, \Delta t\r) \pm \Delta t \Phi \l(t_i,y(t_i),y_\ipo,\Delta t\r).
\end{align*}
Using \eqref{II.3} and \eqref{II.4} we obtain
\begin{align*}
|e(t_\ipo,\Delta t)|
& \leq |e(t_i,\Delta t)| + \Delta t L|y(t_\ipo)-y_\ipo| + \Delta t L|y(t_i)-y_i| - \eta(t_i,\Delta t) \\
& = |e(t_i,\Delta t)| + \Delta t L|e(t_{i+1},\Delta t)| + \Delta t L| e(t_i,\Delta t) | + |\eta(t_i,\Delta t)|.
\end{align*}
This gives
\begin{align*}
(1-\Delta t L) |e(t_\ipo,\Delta t)| & \leq (1+\Delta t L) |e(t_i,\Delta t)| + |\eta(t_i,\Delta t)|,
\end{align*}
so
\begin{align*}
|e(t_\ipo,\Delta t)| & \leq \frac{(1+\Delta t L)}{(1-\Delta t L)} |e(t_i,\Delta t)| + \einsdurch{(1-\Delta t L)}|\eta(t_i,\Delta t)| \\
& \leq \frac{(1+\Delta t L)}{(1-\Delta t L)} |e(t_i,\Delta t)| + \einsdurch{(1-\Delta t L)}\eta(\Delta t).
\end{align*}
By setting
\begin{align*}
\varrho_i & := \frac{(1+\Delta t L)}{(1-\Delta t L)}-1 = \frac{2\Delta t L}{1-\Delta t L} \geq 0 \\
z_i & := |e(t_i,\Delta t)| \geq 0 \\
\eta_i & := \einsdurch{(1-\Delta t L)}\eta(\Delta t) \geq 0
\end{align*}
and applying Lemma II.5 we obtain
\begin{align*}
|e(t_\ipo,\Delta t)|
& = z_\ipo \\
& \leq \l(z_0 + \sum_{k=0}^{i}\eta_k\r)e^{\sum_{k=0}^{i }\varrho_k} \\
& = \l(|e(t_0,\Delta t)| + \sum_{k=0}^{i }\einsdurch{1-\Delta t L}\eta(\Delta t)\r) e^{\sum_{k=0}^{i }\frac{2\Delta t L }{1-\Delta t L}}. \tag{$\ast$}
\end{align*}
Observe that the two sums can be rewritten as
\begin{align*}
\sum_{k=0}^{i }\einsdurch{1-\Delta t L }\eta(\Delta t) = \frac{i+1}{1+\Delta t L }\eta(\Delta t) = \frac{t_{i+1}-t_0}{1+\Delta t L} \cdot \frac{\eta(\Delta t)}{\Delta t}
\end{align*}
and
\begin{align*}
\sum_{k=0}^{i }\frac{2\Delta t L}{1-\Delta t L } = (t_\ipo - t_0) \frac{2L}{1-\Delta t L}.
\end{align*}
Inserting this into ($\ast$) gives the result.
\end{proof}

\begin{theorem}
If a one-step method with Lipschitz conditions \eqref{II.3} and \eqref{II.4} is consistent of order $p\in\N$ for an ODE \eqref{II.1} and if the initial value $y_0$ meets
\[
y_0 = \hat y_0 + O(\Delta t^p),
\]
then the method is convergent of order $p$ with respect to \eqref{II.1} and \eqref{II.2}.
\end{theorem}

\begin{remark}\hfill
\begin{itemize}
     \item The error grows exponentially in time.
     \item If in the underlying ODE the Lipschitz-constant $\hat L$ given by
     \[
     \l|f(t,y_1(t))-f(t,y_2(t))\r|\leq \hat L |y_1-y_2|
     \]
     is large, then $L$ from \eqref{II.4} and \eqref{II.5} will also be large.
     \item If the initial condition of the explicit Euler-method meets
     \[
     y_0 = \hat y + O(\Delta t)
     \]
     then it is convergent of first order with respect to the ODE
     \[
     y'(t) = f(t,y(t)) , y(t_0) = \hat y_0.
     \]
\end{itemize}
\end{remark}

\subsubsection*{Runge-Kutta Methods}
We already know that
\[
y(t_\ipo)-y(t_i) \approx \Delta t f(t,y(t)).
\]
Asking whether a better approximation leads to better convergence leads to the Runge-Kutta methods. Trying with the midpoint rule gives
\begin{align*}
\int_{t_i }^{t_\ipo}f(t,y(t))~dt \approx f \l(t_i + \tfrac{\Delta t }{2},~y\l(t_i + \tfrac{\Delta t}{2}\r)\r)\cdot \Delta t.
\end{align*}
But we have not evaluated $y(t_i,\tfrac{\Delta t }{2})$, what is a problem. The idea is using the explicit Euler-method to approximate $y(t_i+\tfrac{\Delta t}{2})$. Define
\[
y_{i+\teinhalb}:= y_i + \frac{\Delta t }{2} f(t_i, y(t_i)).
\]
By plugging in we obtain
\[
y_{i+1} = y_i + \Delta t \cdot f \l(t_i+\tfrac{\Delta t}{2}, y_i + \tfrac{\Delta t}{2} f(t_i,y(t_i))\r)
\]
which is also called the explicit midpoint rule. \newpage

\subsubsection*{\underline{Excursion to quadrature}}
\definecolor{shadethmcolor}{HTML}{F6FDFF}
We know, that
\[
I: C \to \R, ~ f\mapsto \int_{a }^{b }f(\tau )~d\tau
\]
is a linear functional from some function space $C$ into the real numbers.

\begin{exdef}
A function $Q_{n+1}\in C([a,b])$ with
\[
Q_{n+1}(f) = \sum_{i=0 }^{n}a_i f(x_i)
\]
with nodes $x_i\in[a,b]$ and weights $a_i\in\R$ is called a quadrature rule. Its quadrature error is the linear functional
\[
R_{n+1}(f)=I(f)-Q_{n+1}(f).
\]
The rule converges, if it holds
\[
\limn Q_{n+1}(f) = I(f).
\]
\end{exdef}

\noindent There exist many quadrature rules. Here we consider \textbf{quadrature by interpolation}. Let us assume, that we know $f$ only at $(n+1)$ points $x_0,\dots, x_n$ an we interpolate $f$ by a polynomial $p$ with degree $n$.

\begin{exdef}[Quadrature by interpolation]
Let $p_n$ be a polynomial of degree $n$ on the interval $[a,b]$. We call $R_{n+1}$ a \textbf{quadrature rule by interpolation} if
\[
R_{n+1}(p_n) = 0,
\]
i.e. if a polynomial of degree $n$ can be integrated exactly.
\end{exdef}
\newpage \definecolor{shadethmcolor}{HTML}{EEEEEE}

\noindent Consider $\xi_j = t_i + c_j \Delta t, ~c_j \in [0,1]$ for $j=1,\dots,s$. Then we have
\[
y(t_\ipo) - y(t_i) = \int_{t_i }^{t_\ipo}f(\tau,y(\tau))~d\tau \approx \Delta t \sum_{j=1}^{s }b_j f(\xi_j,y(\xi_j)).
\]
From the quadrature by interpolation we know
\[
\sum_{j=1 }^{s }b_j = 1.
\]
But since we do not know the values $y(\xi_j)$, we have to think about how to get these values.

Applying the fundamental theorem gives
\begin{align*}
y(\xi_j)-y(t_i) &= \int_{t_i }^{t_i + c_j \Delta t}f(t,y(t))~dt\\
& \approx c_j \Delta t \cdot \sum_{\nu=1}^{s} \tilde{a}_{j\nu} f(\xi_\nu,y(\xi_\nu)).
\end{align*}
This seems strange, since we are using the same $\xi_i$. Setting $a_{j\nu}:=c_j \tilde{a}_{j\nu}$ we obtain
\[
k_i = y_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu}f(\xi_\nu,y(\xi_\nu))
\]
as an approximation of $y(\xi_i)$, where $i=1,\dots,s$. These methods are called the \textit{Runge-Kutta methods}.

\begin{definition}[Runge-Kutta-Method/RKM]
For $b_j,c_j,a_{j\nu}\in\R,~j=1,\dots,s$ we denote
\begin{align*}
k_j=y_i + \Delta t \sum_{\nu=1}^{s}a_{j\nu}f(\xi_\nu,k_\nu)
\end{align*}
for $j=1,\dots,s$ and
\begin{align*}
y_{i+1}=y_i + \Delta t \sum_{j=1}^{s} b_j f(\xi_j,k_j)
\end{align*}
with $\xi_j=t_i+c_j\Delta t$ as an \textbf{$s$-step Runge-Kutta method}. We call $c_j$ and $b_j$ \textbf{weights}.
\end{definition}

\begin{remark}\hfill
\begin{itemize}
     \item A Runge-Kutta methods is defined by the parameters $a_{j\nu},b_j,c_j\in\R$.
     \item The \textit{Butcher table} or \textit{Array} of a Runge-Kutta method can be denoted as
     \[
     \begin{matrix}[c|cccc]
     c_1 & a_{11} & \dots & a_{1s} \\ \vdots & \vdots & &\vdots \\ c_s & a_{s1} & \cdots & a_{ss} \\ \hline & b_1 & \cdots & b_s
     \end{matrix}~.
     \]
\end{itemize}
\end{remark}

\begin{example}\hfill
\begin{enumerate}[(1)]
\item For the explicit Euler method, the Butcher table is given by $~\begin{matrix}[c|c]0&0\\ \hline &1\end{matrix}~$, which is equivalent to $k_1=y_i$ and $y_{i+1}=y_i+\Delta t f(t_i,k_1)$.
\item For the implicit Euler method, the Butcher table is given by $~\begin{matrix}[c|c]1&1\\ \hline &1\end{matrix}~$, which is equivalent to $y_{i+1}=k_1=y_i+\Delta t f(t_i+\Delta t,k_1)$.
\item For the explicit midpoint rule, the Butcher tabel is given by
\[
\begin{matrix}[c|cc]
0 & 0 & 0 \\ \teinhalb & \teinhalb & 0 \\ \hline & 0 &1
\end{matrix}
\]
which is equivalent to
\begin{align*}
k_1 &= y_i\\
k_2 &= y_i + \frac{\Delta t }{2 } f(t_1,k_1) \\
y_{i+1} &= y_i + \Delta t f\l(t_i+\tfrac{\Delta t}{2}, k_2\r).
\end{align*}
\end{enumerate}
\end{example}

\noindent The Runge-Kutta method can also be seen from a predictor-corrector method-point of view. For that, consider the trapezoid rule given by
\[
y(t_{i+1})-y(t_i) \approx \frac{\Delta t }{2 }\Big(f(t_i,y(t_i))+f(t_\ipo,y(t_\ipo))\Big).
\]
We obtain
\[
y_\ipo = y_i + \frac{\Delta t }{2} \Big(f(t_i,y_i)+f(t_\ipo,y_\ipo)\Big).
\]
Since we do not know $y_\ipo$, we approximate $f(t_\ipo,y_\ipo)$ by $f(t_\ipo,k_2)$, where
\[
k_2 = y_i + \Delta t f(t_i,y_i)
\]
is derived from the explicit Euler method. We thus have
\begin{align*}
k_1 & = y_i \\
k_2 & = y_i + \Delta t f(t_i,y_i) \\
y_\ipo & = y_i + \frac{\Delta t }{2} \Big(f(t_1,k_1)+f(t_\ipo,k_2)\Big).
\end{align*}
In this case, $k_2$ is called a predictor and since we have
\[
y_{i+1} = k_2 + \frac{\Delta t }{2 }\Big(f(t_i,k_2)-f(t_i,k_1)\Big),
\]
the term
\[
\frac{\Delta t }{2 }\Big(f(t_i,k_2)-f(t_i,k_1)\Big)
\]
is called the corrector. \newline

\noindent Instead of computing $k_j$ at $y(\xi_j)$ we can use the slopes or gradients
\[
r_j = f(t_i+c_j\Delta t, k_j).
\]
Within the predictor-corrector method we have
\begin{align*}
r_1 & = f(t_i,y_i) \\
r_2 & = f\l(t_i + \Delta t, y_i + \Delta t r_1\r)\\
y_\ipo & = y_i + \frac{\Delta t }{2 }(r_1+r_2),
\end{align*}
where $r_1+r_2$ can be seen as an intermediate slope. By setting $r_j = f(t_i+c_j\Delta t, k_j)$ we obtain a Runge-Kutta method
\begin{align*}
r_j & = f(t_i + c_j \Delta t, k_j) \\
& = f\l(t_i + c_j \Delta t, y_i+\Delta t\cdot \sum_{\nu=1}^{s }a_{j\nu}f(\xi_\nu,k_\nu)\r) \\
& = f \l(t_i+c_j\Delta t, y_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu} r_\nu\r).
\end{align*}
By summing up, we can write
\begin{align*}
y_\ipo = y_i + \Delta t \sum_{j=1 }^{s }b_j r_j.
\end{align*}

\noindent When computing for example $r_3$, it can happen that we end up with the form
\[
r_3 = f \l(\dots, \sum_{\nu=1}^{s}a_{j\nu}r_\nu\r),
\]
where $r_3$ depends on $r_3$. Analogue to the one-step methods, in this case we call the Runge-Kutta method implicit. Let us assume, that $A=[a_{j\nu}]\in\R^{s,s}$ is a strict lower triangle matrix, i.e. $a_{j\nu}=0$ for $\nu\geq j$. Then we obtain
\begin{align*}
r_j = f \l(t_i + c_\nu \Delta t, y_i + \Delta t\cdot \sum_{\nu=1}^{j-1}a_{j\nu} r_\nu\r)
\end{align*}
for $j=1,\dots,s$ and hence we have an \textbf{explicit Runge-Kutta method}. If we don't have such a matrix $A$, we have an \textbf{implicit Runge-Kutta method}.

Note, that unlike to one-step methods, in Runge-Kutta methods explicit and implicit only refers to the intermediate steps between $t_i$ and $t_\ipo$. \newline

\noindent Let us assume, that we have a full matrix A and $f:[a,b]\times \R^m \to \R^m$. Then
\begin{align*}
r_1 & = f \l(t_i + c_1 \Delta t, y_i + \Delta t \cdot \sum_{\nu=1}^{s }a_{1\nu} r_\nu\r) \\
&\cdots \tag{II.$\ast$} \\
r_s & = f \l(t_i + c_s\Delta t, y_i + \Delta t \cdot \sum_{\nu=1}^{s }a_{s\nu} r_\nu \r)
\end{align*}
is a system of dimension $s-m$ for computing the gradients $r_j\in\R^m$. It might be linear or non-linear, depending on the underlying system.

\begin{example}[Classical Runge-Kutta method]
Let the Butcher table be given by
\[
\begin{matrix}[c|cccc]
0 & 0 & 0 & 0 & 0 \\
\teinhalb & \teinhalb & 0 & 0 & 0 \\
\teinhalb & 0 & \teinhalb & 0 & 0 \\
1 & 0 & 0 & 1 & 0 \\ \hline
& \tfrac{1}{6} & \tfrac{1}{3} & \tfrac{1}{3} & \tfrac{1}{6}
\end{matrix}~.
\]
This Butcher table gives an explicit Runge-Kutte method, since $A$ is a strict lower triangular matrix. We thus have
\[
y_\ipo = y_i + \Delta t \Big( \tfrac{1}{6} r_1 + \tfrac{1}{3} r_2 + \tfrac{1}{3} r_3 + \tfrac{1}{6} r_4 \Big)
\]
and further
\begin{align*}
r_1 & = f(t_i+ 0\cdot\Delta t, y_i + \Delta t\cdot 0) = f(t_i,y_i) \\
r_2 & = f \l(t_i + \frac{\Delta t}{2}, y_i + \frac{\Delta t }{2} r_1\r) \\
r_3 & = f \l(t_i + \frac{\Delta t }{2}, y_i + \frac{\Delta t }{2 }r_2\r) \\
r_4 & = f(t_1 + \Delta t, y_1 + \Delta t r_3).
\end{align*}
{\rot \texttt{The drawing is missing.}} A Runge-Kutta method hence allows us to approximate $f(t,y(t))$ on intermediate steps, i.e. estimate $y(t)$ on intermediate steps. %This can be very useful, since in many cases the function $y$ is hard or expensive to evaluate.
\end{example}


\noindent Reconsider the (maybe non-linear) system (II.$\ast$) for computing the gradients $r_j$. The following theorem shows that the system can be solved if $f$ satisfies certain conditions. In particular, we can, in some sense, buy the solvability of the system (II.$\ast$) by choosing smaller time steps.

\begin{theorem}
Let the mapping $f:[a,b] \times \R^m \to \R^m$ be continuous so that it holds
\[
\norm{ f(t,\tilde y) - f(t,y) }_\infty \leq L \norm{\tilde y - y}_\infty
\]
for all $t\in [a,b]$, where $L>0$ is a Lipschitz constant. Consider the Runge-Kutta method $(A,b,c)$ with $\Delta t < \tfrac{1}{L\norm{A}_\infty}$. Then for any $j=1,\dots,s$, the iteration given by
\[
r_j^{(\ell+1)} = f \l(t_i + c_i \Delta t, ~ y_i + \Delta t \sum_{\nu=1 }^{s }a_{j\nu} r_\nu^{(\ell)}\r)
\]
converges for any arbitrary initialization $r_1^{(0)}, \dots, r_s^{(0)}$ for $\ell \to \infty$ to the unique solution of the system
\[
r_j= f\l(t_i+c_j\Delta t,~y_i+\Delta t \sum_{\nu=1 }^{s }a_{j\nu} r_\nu\r).
\]
\end{theorem}
\begin{proof}
We set
\[
R:=\begin{bmatrix} r_1 \\ \vdots \\ r_s \end{bmatrix} \und F:= \begin{bmatrix} F_1 \\ \vdots \\ F_2 \end{bmatrix} : \R^s \to \R^s
\]
with
\[
F_j(R) = f \l(t_i + c_j\Delta t,~ y_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu}r_\nu \r).
\]
We thus have
\begin{align*}
\norm{F(R)-F(\tilde R)}_\infty
& \leq L \cdot \norm{
\begin{bmatrix} \Delta t \sum_{\nu=1}^{s}a_{1\nu}(r_\nu - \tilde r_\nu) \\
      \vdots \\
      \Delta t \sum_{\nu=1}^{s}a_{s\nu}(r_\nu - \tilde r_\nu)
\end{bmatrix}}_\infty
\end{align*}
with
\[
\norm{
\begin{bmatrix} \Delta t \sum_{\nu=1}^{s}a_{1\nu}(r_\nu - \tilde r_\nu) \\
      \vdots \\
      \Delta t \sum_{\nu=1}^{s}a_{s\nu}(r_\nu - \tilde r_\nu)
\end{bmatrix}}_\infty \leq \norm{
\begin{bmatrix} \Delta t \sum_{\nu=1}^{s}a_{1\nu} \\
      \vdots \\
      \Delta t \sum_{\nu=1}^{s}a_{s\nu}
\end{bmatrix}}_\infty  \cdot \norm{R - \tilde R}_\infty.
\]
We obtain
\[
\norm{F(R)-F(\tilde R)} \leq L \cdot \Delta t \cdot \underbrace{\l(\max_{i=1,\dots,s} \sum_{\nu=1}^{s}|a_{i\nu}| \r)}_{= \norm{A}_\infty} \norm{R-\tilde R}_\infty
\]
which is a contraction in the Banachspace $(\R^s, \norm{\cdot}_\infty)$ if $L\cdot \Delta t \norm{A}_\infty
 < 1$. The Banach-fixpoint theorem then implies that there exists an $R\in\R^s$ as the fixed point of this iteration. Furthermore, $(R^{(\ell)})_\ell$ with $R^{(\ell+1)} = F(R^{(\ell)})$ converges towards R.
\end{proof}

\begin{remark}
After Definition II.2 (consistency) we saw, that the minimum requirement for a consistent method is
\[
\lim_{\Delta t \to 0} \Phi(t,y(t),y(t+\Delta t),\Delta t) = f(t,y). \labtag{II.6}
\]
\end{remark}

\noindent We now are interested in the consistency of arbitrary Runge-Kutta methods. Reconsider
\begin{align*}
r_j & = f\l(t_i + c_j \Delta t, ~ y(t_i) + \Delta t \sum_{\nu=1}^{s }a_{j\nu }r_\nu\r) \\
& = f(t_i, y(t_i)) + O(\Delta t),
\end{align*}
since
\begin{align*}
y_\ipo & = y_i + \Delta t \cdot \Phi(t_i, y(t_i), y(t_i+\Delta t),\Delta t) \\
& = y_i + \Delta t \cdot \sum_{j=1 }^{s} b_j r_j.
\end{align*}
We further have that
\[
\Phi(t_i, y(t_i), y(t_i+\Delta t), \Delta t) = \sum_{j=1}^{s }b_j f(t,y(t_i))+ O(\Delta t)
\]
and hence
\[
\lim_{\Delta t \to 0} \Phi(t_i, y(t_i), y(t_i+\Delta t), \Delta t)  = f(t_i,y(t_i)) \gdw \sum_{j=1}^{s }b_j = 1.
\]
In summary, we have proven the following result.

\begin{lemma}
A Runge-Kutta method given by $(A,b,c)$ is consistent in the sense of \eqref{II.6} if and only if it holds
\[
\sum_{j=1}^{s }b_j = 1.
\]
\end{lemma}

This makes hope for determining the order of consistency of a Runge-Kutta method by looking at its Butcher table. The next theorem gives us precise conditions to that.

\begin{theorem}
For a Runge-Kutta method $(A,b,c)$ it holds
\begin{enumerate}[(a)]
\item The method has order of consistency $p\geq1$ if
\[
\sum_{j=1}^{s }b_j = 1 \labtag{II.7}
\]
holds for all $j=1,\dots,s$.
\item The method has order of consistency $p\geq2$ if it holds \eqref{II.7} and
\[
\sum_{j=1}^{s }b_j c_j = \einhalb \labtag{II.8}
\]
holds.
\item The method has order of consistency $p\geq3$ if \eqref{II.7}, \eqref{II.8} and
\[
\sum_{j=1}^{s }b_jc_j^2 = \einsdurch{3} \und \sum_{j=1}^{s} b_j \sum_{\nu=1}^{s }a_{j\nu} c_\nu = \einsdurch{6}
\]
hold.
\end{enumerate}
\end{theorem}
\begin{proof}
The proof can be found in \textit{Deuflhard/Bornemann: Numerische Mathematik II, Chapter 4} and is not given in this lecture.
\end{proof}

\begin{remark}\hfill
\begin{itemize}
     \item In the case of an explicit Runge Kutta method the conditions in theorem II.11 are equivalent conditions to the consistency.
     \item The proof of the theorem relies on Taylor expansion.
     \item In principle, this technique is also working for higher orders. However, for higher orders we get more conditions on the coefficients $(A,b,c)$ of the Runge-Kutta method. The number of conditions $N_p$ increases rapidly when the order increases, since for $p=15$ there are already $141 083$ conditions needed.
     \item It seems to be, that the number of stages (steps) of an explicit Runge-Kutta method gives the order. But it holds

     \centering{\begin{tabular}{c|c|c|c|c||c|c|c|c|c}
          order $p$ & $1$ & 2 & 3 & 4 & 5 & 6 & $\cdots$\!& 12 & 14 \\ \hline
          $s$ (stages) & 1 & 2 & 3 & 4 & 6 & 7 &  $\cdots$\!& 25 & 35
     \end{tabular}.}
\end{itemize}
\end{remark}

\begin{theorem}
For an explicit Runge-Kutta method with $s$ stages and order of consistency $p$ it holds
\[
p \leq s.
\]
\end{theorem}
\begin{proof}
If we can find one IVP such that $p>s$ leads to a contradiction, the claim is proven. Consider
\[
y'(t) = y(t), ~ y(0) = 1.
\]
\end{proof}

\noindent In chapter I we saw that any non-autonomous system
\begin{align*}
y'(t) & = f(t,y(t)) \\
y(t_0) & = y_0
\end{align*}
can be transformed to an autonomous system
\begin{align*}
z'(t) & = \hat f(z(t)) \labtag{II.9}\\
z'(t_0) & = z_0 = \begin{bmatrix} y_0 \\ t_0 \end{bmatrix}
\end{align*}
by the transformation
\begin{align*}
z & =\begin{bmatrix} y \\ s \end{bmatrix} \\
\hat f(z) &= \begin{bmatrix} f(s,y) \\ 1 \end{bmatrix}.
\end{align*}
We want to know if the result of an explicit Runge-Kutta method stays the same when we apply the Runge-Kutta method to the autonomous system.

By applying the explicit Runge-Kutta method to the autonomous sytem \eqref{II.9}, we obtain
\begin{align*}
z_{i+1} & = z_i + \Delta t \sum_{j=1}^{s} b_j \hat r_j,
\end{align*}
or more precisely
\begin{align*}
\begin{bmatrix} y_{i+1} \\ t_{i+1} \end{bmatrix}  & = \begin{bmatrix} y_i \\ t_i  \end{bmatrix}  + \Delta t \sum_{j=1}^{s }b_j \begin{bmatrix} r_j \\ 1 \end{bmatrix}.
\end{align*}
We will first just consider the lower part of the last equation, i.e.
\[
t_{i+1} = t_i + \Delta t \sum_{j=1}^{s }b_j = t_i + \Delta t,
\]
so it is necessery that
\[
\sum_{j=1}^{s }b_j = 1.
\]
Further, we have that
\begin{align*}
\hat r_j & = \begin{bmatrix} r_j \\ 1 \end{bmatrix} = \hat f \l(z_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu} \hat r_\nu\r) \\
& = \hat f \l(\begin{bmatrix} y_i \\ t_i  \end{bmatrix} + \Delta t \sum_{\nu=1}^{s }a_{j\nu } \begin{bmatrix} r_j \\ 1 \end{bmatrix} \r),
\end{align*}
such that\footnote{Remember that in the non-autonomous case the structure of $r_j$ is
\[
r_j = f \l(t_i + c_j \Delta t, ~ y_i + \Delta t \sum_{\nu=1 }^{s }a_{j\nu} r_\nu\r).
\]
}
\begin{align*}
r_j = f \l(t_i+\Delta t \sum_{\nu=1}^{s }a_{j\nu}\cdot 1, ~ y_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu }r_\nu\r),
\end{align*}
so our second condition is
\[
\sum_{\nu=1 }^{s }a_{j\nu} = c_j.
\]
We summarize our result in the following lemma.

\begin{lemma}
The explicit Runge-Kutta method $(A,b,c)$ is invariant to making it autonomous if and only if it holds
\begin{enumerate}[(i)]
\item
\[
\sum_{j=1 }^{s }b_j = 1
\]
\item
\[
\sum_{\nu=1}^{s }a_{j\nu} =c_j
\]
for all $j=1,\dots,s$.
\end{enumerate}
\end{lemma}

We are heading towards the implicit Runge-Kutta methods and try to motivate them.

Since in the implicit case the matrix $A$ has more non-zero entries, we have more conditions to our method. We seek to find out if it is possible to use the additional coefficients of an implicit Runge-Kutta method to obtain a higher order for a given number $s$ of stages. In the explicit case, theorem II.12 shows that for the order of consistency $p$ it holds $p\leq s$, so we want to find out if it is possible that $p>s$ for a implicit Runge-Kutta method.

Our idea is to use collocation method, i.e. choose (simple) functions from some function space (e.g. polynomials $\mathbb{P}$) and a set of collocation points (pairwise different points) and set the free coefficients of the (simple) function such that the problem function (e.g. $f$) holds in the collocation points. Consider the ODE
\[
y'(t) = f(t,y(t)), ~y(t_0)=y_0.
\]
We assume that the numerical solution has been carried up to the point $(t_i,y_i)$. We now seek for a recipe to advance it to $(t_{i+1},y_{i+1})$, where $t_{i+1}=t_i+\Delta t$. To do so, we choose $s$ collocation points $c_1,\dots,c_s\in[0,1]$ and then seek for a $s$-th degree polynomial $u_s\in\mathbb{P}_s$, such that it holds
\begin{align*}
u_s(t_i) & = y_i \labtag{II.10} \\
u_s'(t_i+c_j\Delta t) & = f\Big(t_i+c_j\Delta t,~u_s(t_i+c_j\Delta t)\Big)
\end{align*}
for $j=1,\dots,s$.

A \textbf{collocation method} consists of finding such a $u_s$ and setting $y_{i+1} = u(t_{i+1})$.

But we do not yet know what the relation to a Runge-Kutta method is.

\begin{lemma}
The collocation method defined by \eqref{II.10} is equivalent to a $s$ stage implicit Runge-Kutta method with the coefficients
\begin{align*}
a_{ji}& =\int_{0}^{c_j }L_i(\tau)~d\tau \labtag{II.11} \\
b_j & = \int_{0}^{1}L_j(\tau )~d\tau,\labtag{II.12}
\end{align*}
where $L_j(\tau)$ is the Lagrangian interpolation polynomial
\[
L_j(\tau) = \prod_{\substack{\ell = 1 \\ \ell \neq j}}^{s} \frac{\tau - c_\ell}{c_j - c_\ell}.
\]
\end{lemma}
\begin{proof}
By the collocation polynomial $u_s(t)$ we define $r_j:= u_s'(t_i+c_j\Delta t)$. By the Lagrange interpolation formula, for any $\tau \in [0,1]$ we have that
\[
u'_s(t_i+\tau\Delta t) = \sum_{\ell=1}^{s }L_\ell(\tau)r_\ell.
\]
Integration gives
\begin{align*}
u_s(t_i+c_j\Delta t) & = u_s(t_i) + \Delta t \cdot \sum_{\ell=1 }^{s }r_\ell \int_{0}^{c_j}L_\ell(\tau)~d\tau
\end{align*}
for all $j=1,\dots,s$. Inserting for \eqref{II.11} gives
\[
u_s(t_i+c_j\Delta t) = u_s(t_i)+\Delta t  \sum_{\ell=1}^{s}r_\ell a_{j\ell}.
\]
Since
\[
r_j = u_s'(t_i+c_j\Delta t)\overset{\text{\eqref{II.10}}}{=} f\Big(t_i+c_j\Delta t,~u_s(t_i+c_j\Delta t)\Big),
\]
we obtain
\[
u_s(t_i+c_j\Delta t) = u_s(t_i) + \Delta t \sum_{\ell=1}^{s }a_{j\ell} f\Big(t_i+c_j\Delta t,~u_s(t_i+c_j\Delta t)\Big),\tag{$\ast$}
\]
or
\[
u_s(t_i+\Delta t) = u_s(t_i) + \Delta t \sum_{j=1}^{s }b_j f\Big(t_i+c_j\Delta t,~u_s(t_i+c_j\Delta t)\Big). \tag{$\ast\ast$}
\]
Since in the collocation method we set $u_s(t_i)= y_i$ and $u_s(t_{i+1})=y_{i+1}$, we have the Runge-Kutta method in ($\ast$) and ($\ast\ast$).
\end{proof}

We do not know yet if every Runge-Kutta method originates in collocation. In general, this is not true, as the following example demonstrates: Consider the two stage Runge-Kutta method with $c_1=0,~c_2=\tfrac{2}{3}$. Computing $L_1,L_2$ and $a_{ji},b_i$ via \eqref{II.11} and \eqref{II.12} will give
\[
\begin{matrix}[c|cc]
     0 & 0 & 0 \\
     \tfrac{2}{3} & \tfrac{1}{3} & \tfrac{1}{3} \\ \hline
     & \tfrac{1}{4} & \tfrac{3}{4}
\end{matrix}~~.
\]
Given that every choice of collocation points corresponds to an unique collocation method (Lagrange interpolation, Numerical Mathematics I) we deduce, that the implicit Runge-Kutta method
\[
\begin{matrix}[c|cc]
     0 &\tfrac{1}{4}& -\tfrac{1}{4} \\
     \tfrac{2}{3} & \tfrac{1}{4} & \tfrac{5}{12} \\ \hline
     & \tfrac{1}{4} & \tfrac{3}{4}
\end{matrix}
\]
with order $p\geq 3$ has no collocation counterpart.

\begin{example}
For $s=1$ the polynomial defined by
\[
u(t) = y_i + (t-t_i)f(t_i+c_1\Delta t,~u(t_1+c_1\Delta t))
\]
gives a collocation method like in \eqref{II.10}. For $c_1 = 0$ it gives the explicit Euler method, for $c_1 = 1$ it gives the implicit Euler method and for $c_1=\teinhalb$ the midpoint rule is given. {\rot \texttt{The drawing is missing.}}
\end{example}

We want to know how the order $p$ of an implicit Runge-Kutta method depends on the choice of the collocation method. To answer this question reconsider \eqref{II.1}, i.e. the ODE
\[
y'(t) = f(t,y(t))
\]
and let $v$ be a "candidate" solution. The defect is defined as
\[
d(t,v) = v'(t) - f(t,v(t)). \labtag{II.13}
\]
When the defect is small, the error $\norm{y(t) - v(t)}$ is small too, because if $d(t,v)=0$ our candidate $v$ is the solution.

Let us consider the linear case
\[
y' = \Lambda y~, y(t_0)=y_0.
\]
Then the defect is given by
\[
D(t) = v'(t) - \Lambda v(t).
\]
Plugging in for our candidate solution $v$ we thus arrive at the linear inhomogeneous ODE
\[
v'(t) = \Lambda v(t) + D(t), ~t\geq t_0, ~v(t_0) \textrm{ given.}
\]
The exact solution of this ODE (theorem I.12) is given by
\[
v(t) = e^{(t-t_0)\Lambda} + \int_{t_0}^{t }e^{(t-\tau)\Lambda} D(\tau)~d\tau
\]
wheras the solution of \eqref{II.13} is given by
\[
y(t) = e^{(t-t_0)\Lambda}y_0 ~ t\geq t_0.
\]
Hence
\begin{align*}
v(t) - y(t)
& = e^{(t-t_0)\Lambda}(v_0-y_0) + \int_{t_0}^{t}e^{(t-\tau)\Lambda}D(\tau)~d(\tau)
\end{align*}
for $t\geq t_0$. Thus, the error can be controlled by the observables $v_0, y_0$ and the defect $D(\tau)$.

\begin{theorem}[Alekseenz-Gröbner-Lemma]
Let $v$ be a smoothly differentiable function that obeys the initial conditions $v(t_0)=y_0$. Then
\[
v(t) - y(t) = \int_{t_0}^{t }\Xi(t,\tau)d(\tau,v)~d\tau,  ~~~t\geq t_0,
\]
where $\Xi$ is the matrix of the partial derivatives of the solution of the ODE
\[
w'=f(t,w),~w(\tau)=v(\tau).
\]
\end{theorem}
\begin{proof}
The proof can be found in \textit{Deuflhard-Bornemann, Numerical Mathematics, Chapter III, Theorem 3.4} and is not given in this lecture.
\end{proof}

\newpage
\subsubsection*{\underline{Excursion II.2} (Quadrature)}
\definecolor{shadethmcolor}{HTML}{F6FDFF}

\begin{exdef}
A quadrature rule $Q_{n+1}[f]$ is exact of order $n\in\N$ if
\[
R_n[p] = I[p] - Q_{n+1}[p] = 0
\]
holds for all $p \in \P_{n+1}$.
\end{exdef}
Recall that the scalar product of to functions is given by
\[
\scal{f}{g} = \int_{a }^{b }f(x)g(x)w(x)~dx,
\]
where $w(x)$ is a weight-function.
\begin{exthm}
Let $r=\P_s$ obey
\[
\scal{r}{p} = 0
\]
for every $p\in \P_s$ and
\[
\scal{r}{t^s} \neq 0.
\]
Let $c_1,\dots,c_s$ be the zeros of the polynomial $r$ and choose $b_1,\dots,b_s$ such that
\[
\sum_{j=1}^{s }b_j c_j^{s-1} = \int_{a }^{b }\tau^{s-1}~d\tau.
\]
Then the quadrature rule has exactness of order $2s$.
\end{exthm}
\newpage

\definecolor{shadethmcolor}{HTML}{EEEEEE}

\begin{theorem}
An implicit Runge method $(A,b,c)$ generated by the collocation method \eqref{II.10} has for $f\in C^p(\R\times\R^d,\R^d)$ the consistency order $2s$ if and only if the quadrature formula defined by the nodes $c$ and weights $b$ meet the requirements of the previous excursion theorem, i.e. if it is of order $2s$.
\end{theorem}
\begin{proof}
We only give a sketch of the proof. Use theorem II.15 and $u(t_i)=y_i$, where $y_{i+1}$ is the solution of the collocation method and $y(t_{i+1})$ is the true solution. We obtain
\[
y_{i+1}-y(t_{i+1}) = \int_{t_i}^{t_{i+1}}\Xi(t_{i+1},\tau)d(\tau,u)~d\tau.
\]
Applying a quadrature rule with $w(t)=1$ and the quadrature points $t_i+c_1\Delta t, \dots, t_i + c_s\Delta t$ leads to
\[
y_\ipo - y(t_\ipo) = \sum_{j=1}^{s } b_j d(t_i+c_j\Delta t, u) + \textrm{error of quadrature from \eqref{II.14}}.
\]
According to \eqref{II.10} it holds
\begin{align*}
d(t_i+c_j\Delta t, u)
& =  u'(t_i+c_j\Delta t) - f(t_i+c_j\Delta t, u(t_i+c_j\Delta t))\\
{\tiny\textrm{\eqref{II.10}}} & = 0.
\end{align*}
Thus, according to excursion theorem the order of quadrature with the weight functions $w(t)=1, ~0\leq t \leq 1$ and $c_1,\dots, c_s$ is $2s$.
\end{proof}

We have seen that for an explicit Runge-Kutta method it is not possible to choose $p>s$. However, when dealing with implicit Runge-Kutta methods, by a suitable choice of $c_i,a_{ji}$ the implicit Runge-Kutta method can be reffered to as a collocation method. The theorem of Alekseenz-Gröbner tells us that
\[
y_\ipo - y(t_\ipo) \leq \textrm{quadrature error with certain exactness.}
\]
So if we choose the $c_j$ such that they refer to an implicit Runge-Kutta method and also have an quadrature error of order $2s$, then we can obtain up to $p\leq 2s$.

\begin{example}[Gauß-Legendre Runge-Kutta method]
$s$-stage, order $2s$ methods, choose $c_1,\dots,c_s\in(0,1)$ as zeros of the Legendre polynomial $p_s$ given by
\[
p_s(t) = \frac{(s!)^2}{(2s)!} \cdot \sum_{k=0}^{s} (-1)^k \binom{s}{k} \binom{s+k}{k}t^k.
\]
For $s=1$ we get
\[
p_1(t) = t- \teinhalb,
\]
so $c_1 = \teinhalb$ and the Butcher table is
\[
\begin{matrix}[c|c]
    \teinhalb & \teinhalb \\\hline & 1
\end{matrix}~~.
\]
For $s=2$ we have
\[
p_1(t)=t^2 -t + \tfrac{1}{6},
\]
so $c_1 = \tfrac{1}{2} - \tfrac{\sqrt{3}}{6}$ and $c_2=\tfrac{1}{2}+\tfrac{\sqrt{3}}{6}$ and by \eqref{II.11} and \eqref{II.12} from theorem II.14 we have the Butcher table given by
\[
\begin{matrix}[c|cc]
\teinhalb - \tfrac{\sqrt{3}}{6} & \tfrac{1}{4} & \tfrac{1}{4}-\tfrac{\sqrt{3}}{6} \\
\teinhalb + \tfrac{\sqrt{3}}{6} & \tfrac{1}{4} - \tfrac{\sqrt{3}}{6} & \tfrac{1}{4} \\\hline
& \tfrac{1}{2} & \tfrac{1}{2}
\end{matrix}~~.
\]
\end{example}

From Lemma II.10 we already know that Butcher trees are a method for finding the conditions on $(A,b,c)$ for reaching a certain order of an explicit Runge-Kutta method.  The order of a tree is defined as the number of nodes in the tree. It grows exponentially in the treesize. {\rot \texttt{The tree drawings are missing.}} We have the polynomial $\Phi$ defined by
\[
\Phi(t) = \sum_{i,j} b_i c_i^2 a_{ij} c_j^2
\]
and $\gamma(t) \in \N\ohnenull$ by $\dots$ We obtain the following table:


{\hfill\begin{tabular}{cccc}
tree & order & $\gamma$ & condition \\ \hline
(a) & 1 & 1 & $\sum_{i} b_i = 1$ \\
(b) & 2 & 2 & $\sum_{i} b_i c_i = \teinhalb = \tfrac{1}{\gamma(t)}$ \\
(c) & 3 & 3 & $\sum_{i} b_i c_i^s = \tfrac{1}{3}$
\end{tabular}.\hfill}

\subsubsection*{Stability (A Stability)}
Looking at the explicit Euler method for the solution of the IVP
\[
y'(t) = -20 y(t), ~ y(0)=1.2
\]
one can observe oscillation for large $t$. This motivates the definition of stability of a numerical scheme.

Consider a vector valued IVP given by
\begin{equation}
    \begin{aligned}
        y'(t) & = f(t,y(t)) \\
        y(0) & = \hat y_0
    \end{aligned} \labtag{II.15}
\end{equation}
We want to know what happens with small perturbations of $u_i$ at some $t_i$. We have
\[
(y+u)'(t) = f(t,(y+u)(t))
\]
with
\[
(y+u)(t_i) = y(t_i) + u_i
\]
for $t\geq t_i$. Using Taylor expansion gives
\begin{align*}
u'(t) & = (y+u)'(t) - y'(t) \\
& = f(t,(y+u)(t)) - f(t,y(t)) \\
& \approx f(t,y(t)) + \partiell{f}{y}(t,y(t)) \cdot u(t) - f(t,y(t)) \\
& = \partiell{f}{y} (t,y(t))\cdot u(t).
\end{align*}
Freezing the matrix $\partiell{f}{y} (t,y(t))$ at $t_i$ we obtain a linear system of ODEs with constant coefficients:
\begin{align*}
u'(t) & = \partiell{f}{y} (t,y(t_i)) \cdot u(t)
\end{align*}
or
\begin{align*}
u'(t) & = A u(t),
\end{align*}
where $A$ is a matrix.

Remember, that for an $A\in\C^{n,n}$ the ODE
\[
y'(t) = A y(t) \tag{$\ast$}
\]
has for $n=1$ the solution
\[
y(t) = y(t_0)e^{at},
\]
where $a=A$. Using this as an Ansatz for ($\ast$) we obtain
\[
y(t) = v e^{\lambda t},~\lambda\in \C, v\in\C^n.
\]
Hence, $y(t)=ve^{\lambda t}$ solves ($\ast$) if and only if $\lambda$ is an eigenvalue of $A$ and $v$ is an eigenvector of $A$, since
\[
\lambda v e^{\lambda t} = A v e^{\lambda t}.
\]

From Analysis II, recall that if $A$ is diagonalizable, then the matrix
\[
Y = (v_1 e^{\lambda_1},\dots,v_n e^{\lambda_n})
\]
is the fundamental matrix of $y'(t) = Ay(t)$. \newline

\noindent Heading back, this means that
\[
u(t) = c_1v_1 e^{\lambda_1 t } + \dots + c_n v_n e^{\lambda_n t},
\]
where $\lambda_i$ are eigenvalues and $v_i$ are eigenvectors. If $\real(\lambda_i)<0$, then it holds
\[
\lim_{t\to\infty} u(t) = 0.
\]
In this case, we consider the IVP as moderate, since small local perturbations fall off.\footnote{Note, that for sakes of simplicity we assumed pairwise different eigenvalues. However, above standing considerations also hold for multiple eigenvalues.}

We thus consider the scalar test problem
\[
\begin{aligned}
    y'(t) & = \lambda y(t) \\
    y(0) & = 1.
\end{aligned}\labtag{II.D}
\]
This problem is called the \textbf{Dahlquist Test-Problem}.

\begin{definition}[A-stability]
We call a numerical method \textbf{A-stable} (absolute stable) if its approximates $y_i$ of \eqref{II.D} for each $\lambda \in \C^-=\l\{\lambda\in\C ~|~ \real(\lambda)<0\r\}$ with arbitrary but fixed $\Delta t$ are contractive, i.e. if
\[
|y_{i+1}| < |y_i| \labtag{II.16}
\]
for all $i\in\N$
\end{definition}

\begin{example}
\hfill
\begin{enumerate}[(a)]
    \item The explicit Euler method was given by $y_{i+1} = y_i + \Delta t f(t_i,y(t_i))$. Applying this on \eqref{II.D} we obtain
    \begin{align*}
    y_{i+1} & = y_i + \Delta t \cdot\lambda\cdot y_i \\
    & = (1+\Delta t\lambda) y_i,
    \end{align*}
    so we have to choose $\lambda \Delta t \in \l\{z\in\C~|~|z+1|<1\r\}$, since otherwise $(1+\Delta t \cdot \lambda)\not< 1$. Hence, the explicit Euler method is \textbf{not stable}.

    \item The implicit Euler method was given by $y_{i+1}=y_i+\Delta t f(t_{i+1},y_{i+1})$. Applying this to \eqref{II.D} yields
    \begin{align*}
    y(t_{i+1}) = y_i + \Delta t \cdot \lambda y_{i+1},
    \end{align*}
    or equivalently
    \[
    y_{i+1} = \einsdurch{(1-\Delta t \lambda)} y_i.
    \]
    Hence, the implicit Euler method is A-stable.

    \item For the Runge-Kutta methods we have
    \[
    y_{i+1} = y_i + \Delta t \sum_{j=1}^{s }b_j f(t_i+c_j\Delta t, k_j).
    \]
    Applying this to \eqref{II.D} yields
    \[
    y_i + \Delta t \lambda \sum_{j=1}^{s }b_j k_j.
    \]
    Remember that $k_j$ was given by
    \[
    k_j = y_i + \Delta t \sum_{\nu=1}^{s }a_{j\nu}f(t_i+c_j\Delta t,k_\nu) \overset{\tiny\textrm{\eqref{II.D}}}{=}y_i + \Delta t \lambda \sum_{\nu=1}^{s }a_{j\nu} b_j.
    \]
    Rewriting with
    \[
    k= \begin{bmatrix} k_1 & \cdots & k_s \end{bmatrix}^T, ~ e=\begin{bmatrix} 1 & \cdots & 1 \end{bmatrix}^T\in\R^s
    \]
    gives $k=y_i e + \Delta t A k$, where $A=[a_{j\nu}]$. This is equivalent to
    \[
    (I_s - \Delta t \lambda A)k = y_i e.
    \]
    By assuming that $\einsdurch{\Delta t k}\not\in \sigma(A)$, where $\sigma(A)$ is the spectrum of $A$, it holds
    \[
    k=(I_s-\Delta t \lambda A)\inv y_i e.
    \]
    With this, we obtain
    \begin{align*}
    y_{i+1}
    & = y_i+\Delta t \lambda b^T k \\
    & = y_i + \Delta t \lambda b^T (I_s - \Delta t \lambda A)\inv y_i e \\
    & = \underbrace{\l(1+\Delta t \lambda b^T \l(I_s - \Delta t \lambda A\r)\inv e \r)}_{=:R(\lambda\Delta t)}y_i.
    \end{align*}
\end{enumerate}
\end{example}

The observation for Runge-Kutta methods motivates the following definition.

\begin{definition}[Stability function of Runge-Kutta methods]
For $\hat\sigma(A)=\l\{\lambda\in\C\ohnenull ~|~ \lambda\inv \in \sigma(A)\r\}$ we call the mapping
\[
R: \C\setminus \hat\sigma(A) \to \C, ~~\xi\mapsto R(\xi)=1+\xi b^T (I_s-\xi A)\inv e
\]
\textbf{stability function of the Runge-Kutta method} defined by $(A,b,c)$.
\end{definition}

\begin{theorem}
A Runge-Kutta method is A-stable if and only if $|R(\xi)|<1$ for $\xi \in \C^-$, where $R$ is its stability function.
\end{theorem}
\begin{proof}
No proof.
\end{proof}

We want to know how exactly the stability function of a Runge-Kutta method looks like, if we consider explicit or implicit Runge-Kutta methods.

\begin{theorem}
The stability function of a Runge-Kutta method given by $(A,b,c)$ is
\begin{enumerate}[(a)]
    \item a polynomial of degree less or equal $s$ in the case of an explicit RKM and
    \item a rational function, which possibly contains the inverse of the eigenvalues of $A$ as poles in the case of an implicit RKM.
\end{enumerate}
\end{theorem}
\begin{proof}
\hfill
\begin{enumerate}[{ad} (a).]
    \item $A\in\R^{s,s}$ is a strictly lower triangular matrix, hence $A^s = 0$. We consider the Neumann series given by
    \[
    (I_s - \xi A)^{-1} = I_s + \xi A + \dots + \xi^{s-1} A^{s-1}
    \]
    such that
    \[
    R(\xi) = 1 + \xi b^T (I_s \xi A)\inv e = 1+b^T(\xi I_s + \dots +\xi^s A^{s-1})e
    \]
    is a polynomial.
    \item We solve the linear system $(I-\xi A)v = e$ by using Cramer's rule and obtain
    \begin{align*}
    v_j & = \frac{\det \Big( (I-\xi A)_1, \cdots\!, e, \cdots\!, (I-\xi A)_s \Big)}{\det (I-\xi A)}\\
    & = \frac{p_j(\xi)}{\det(I-\xi A)},
    \end{align*}
    with $p_j\in\P_{s-1}$. Thus, we have
    \[
    R(\xi) = 1 + \frac{\sum_{j=1}^{s}b_j p_j(\xi)}{\det(I-\xi A)}\cdot \xi,
    \]
    what ends the proof.
\end{enumerate}
\end{proof}

\begin{theorem}
There exists no A-stable explicit Runge-Kutta method.
\end{theorem}
\begin{proof}
In the explicit case we have $\lim_{|\xi|\to\infty} R(\xi) = \infty$.
\end{proof}

We have thus seen that for explicit RKM the stepsize can not be chosen arbitray. So we seek for $\Delta t$ such that $|R(\lambda \Delta t)|<1$. This leads us to the next definition.

\begin{definition}[Domain of stability]
The set $S:= \l\{\xi \in \C ~|~ |R(\xi)|<1 \r\}$ is called the \textbf{domain of stability} of the RKM belonging to $R$.
\end{definition}

There might be some scenarios where the usage of a constant $\Delta t$ might be suboptimal, since the solution can behave very differently at different times. Thus, in the following, we are not going to assume that $\Delta t$ is constant. Not only the discretization error, but also the \textbf{rounding error} is important. When, by rounding, we obtain $\tilde{y_i}$ instead of $y_i$, the rounding error is given by
\[
|y_i - \tilde y_i| = O \l(\frac{\varepsilon}{\Delta t}\r).
\]
The discretization error is given by
\[
|y(t_i) - y_i| = O(\Delta t^p).
\]
Since the rounding error rises as $\Delta t \to 0$ and the discretization error rises as $\Delta t \to \infty$, we want to know if there exists an optimal stepsize $t_{\textrm{opt}}$.

But how is it possible to find a suitable step size, when the exact solution is unknown? Let $y_i$ be an approximation of $O(\Delta t^p)$ and $\tilde y_i$ be an approximation of $O(\Delta t^q)$ where $q>p$. Then we have\footnote{Note that $y(t_i)$ is unknown.}
\[
\norm{y(t_i) - y_i} \leq \norm{y(t_i)-\tilde y_i} + \norm{\tilde y_i + y_i} = O(\Delta t^q) + \norm{\tilde y_i - y_i}
\]
and hence $O(\Delta t^q) + \norm{\tilde y_i - y_i}$ is an estimator for the exact solution of order $q$.

Consider an equidistant grid of points and a one-step method
\[
\begin{aligned}
y_{\Delta t}(t+\Delta t) & = y_{\Delta t}(t) + \Delta t \Phi(t,y_{\Delta t}(t),\Delta t) \\
y_{\Delta t}(t_0) &= y_0
\end{aligned}\labtag{II.16}
\]
with order of consistency $p\in\N$ and local error (Def II.2)
\[
\eta(t+\Delta t) = y(t) - \Delta t \Phi\big( t,y(t),\Delta t \big) - y(t+\Delta t) = O(\Delta t^{p+1}).\labtag{II.17}
\]
If $\Phi$ is sufficiently often differentiable, we can expand the error and obtain
\[
\eta(t,\Delta t) = \sum_{i=p+1}^{N+1}g_i(t) \Delta t^i + O(\Delta t^{N+2}). \labtag{II.18}
\]
The next theorem shows that a similar expansion is possible for the global error $e(t,\Delta t)$.

\begin{theorem}[Gregg]
Let $f(t,y)$ and $\Phi$ be sufficiently often continuously differentiable on $S=\l\{(t,y) ~|~ t_0 \leq t \leq t_{\textrm{end}},~y\in\R^d\r\}$ and let us assume that the local error $\eta$ has the expansion \eqref{II.18}. Then the global error $e(t,\Delta t)$ after $n$-steps with stepsize $\Delta t$ has in $t^*=t_0 + n\Delta t$ has an asymptotic expansion of the form
\[
e(t,\Delta t) = e_p(t^*)\Delta t^p + e_{p+1}(t^*)\Delta t^{p+1} + \dots + e_N(t^*)\Delta t^N + E_{N+1} (t^*,\Delta t) \Delta t^{N+1},
\]
where the remainder $E_{N+1}(t^*,\Delta t)$ is bounded for $0<\Delta t < \Delta t_0$.
\end{theorem}
\begin{proof}
Observe that
\[
y(t+\Delta t) = y(t) + \int_{t }^{t+\Delta t}f(\tau,y(\tau ))~d\tau \rightsquigarrow y_{\Delta t} (t+\Delta t) = y(t) + \Delta t\Phi(t,y(t),\Delta t),
\]
where
\[
\Phi(t,y(t),\Delta t) = \einsdurch{\Delta t} \int_{t }^{t+ \Delta t }f(\tau, y(\tau))~d\tau - \frac{\eta(t,\Delta t)}{\Delta t}.
\]
We now set
\[
\Phi_y(t,y(t),\Delta t) = f_y(t,y)d\tau + O(\Delta t) \labtag{II.19}.
\]
We construct a method of order $p+1$ such that
\[
y(t) - y_{\Delta t} (t) = e_p(t) \Delta t^p + O(\Delta t^{p+1}).
\]
We set
\[
\bar y_{\Delta t}(t) = e_p(t)\Delta t^p + y_{\Delta t}(t) \labtag{II.20}.
\]
For the "new" one step method we thus have
\[
\begin{aligned}
\bar y_{\Delta t}(t+\Delta t) &= \bar y_{\Delta t} (t) + \Delta t \bar\Phi \l(t, \bar y_{\Delta t}(t), \Delta t\r) \\ \bar y_{\Delta t} (t_0) & = y_0.
\end{aligned} \labtag{II.21}
\]
We then have
\begin{align*}
\bar y_{\Delta t}(t+\Delta t) &\overset{\textrm{\eqref{II.20}}}{=} e_p(t+\Delta t) \Delta t^p + y_{\Delta t}(t + \Delta t) + O(\Delta t^{p+1}) \\
& \overset{\textrm{\eqref{II.16}}}{=} e_p(t+\Delta t) \Delta t^p + y_{\Delta t}(t) + \Delta t \Phi(t,y(t),\Delta t)  + O(\Delta t^{p+1})\\
& ~~~\!=~~ e_p(t+\Delta t)\Delta t^p \overset{\textrm{\eqref{II.20}}}{+} \bar y_{\Delta t}(t) - e_p(t) \Delta t^p \\ & ~~~~~~+ \Delta t \Phi\Big(t,\bar y_{\Delta t}(t) - e_p(t)\Delta t,\Delta t\Big)  + O(\Delta t^{p+1})\\
& ~~~\!=~~ \bar y_{\Delta t}(t) + \Delta t \Phi\Big(t,\bar y_{\Delta t}(t) - e_p(t)\Delta t, \Delta t\Big) \\ & ~~~~~~+ \Big( e_p(t+\Delta t) - e_p(t) \Big)\cdot \Delta t^{p} + O(\Delta t^{p+1}).
\end{align*}
Comparing \eqref{II.21} with this result gives
\[
\Phi\Big(t,y_{\Delta t}(t) - e_p(t)\Delta t, \Delta t\Big) + \Big( e_p(t+\Delta t) - e_p(t) \Big)\Delta t^{p-1} = \bar \Phi ( t, y_{\Delta t}(t),\Delta t).
\]
For $\bar \eta$ it holds
\begin{align*}
\bar \eta (t,\Delta t)
& = y(t) + \Delta t \bar \Phi(t,y(t),\Delta t) - y(t+\Delta t) \\
& = y(t) + \Delta t \Phi\Big(t,y(t)-e_p(t)\Delta t^p, \Delta t\Big) \\
& ~~~+ \Big(e_p(t+\Delta t) - e_p(t)\Big) \Delta t^p - y(t+\Delta t) + O(\Delta t^{p+1}) \\
& = y(t) + \Delta t \Phi\Big(t,y(t)-e_p(t)\Delta t^p, \Delta t\Big) +  e_p' \Delta t^{p+1} - y(t+\Delta t) + O(\Delta t^{p+1}) \\
& = y(t) + \Delta t \Phi\big(t,y(t),\Delta t\big)-f_y(t,y)\cdot e_p(t)\Delta t^{p+1} \\
& ~~~+  e_p'(t) \Delta t^{p+1} - y(t+\Delta t) + O(\Delta t^{p+1}) \\
& = \Delta t^{p+1} \Big( g_{p+1}(t) + e_p'(t) - f_y(t,y(t)) e_p(t) \Big) + O(\Delta t^{p+2})
\end{align*}
and thus the method \eqref{II.21} has order of consistency $p+1$, if $e_p(t)$ solves the IVP
\begin{align*}
e_p(t) & = f_y(t,y(t)) e_p(t) + g_{p+1}(t) \\
e_p(t_0) & = 0.
\end{align*}
For the global error of \eqref{II.21} it holds (according to theorem II.7)
\[
y(t^*) - \bar y_{\Delta t}(t^*) = E_{p+1}(t,\Delta t) \Delta t^{p+1}
\]
and with \eqref{II.20} it holds
\[
e(t^*,\Delta t) = e_p(t^*)\Delta t^p + E_{p+1} (t^*,\Delta t)\Delta t^{p+1}.
\]
For the next step $e_{p+1}(t)$ we start the proof at \eqref{II.21}.
\end{proof}

With theorem II.23 in our hands, we can
\begin{enumerate}[(a)]
    \item construct a method of higher order and
    \item construct an error estimator for the global error.
\end{enumerate}
\begin{enumerate}[{ad} (a):]
    \item Let
        \begin{align*}
        y_{i+1}^{\Delta t}&, ~\textrm{stepsize}~ \Delta t, ~ 1 \textrm{ step}\\
        y_{i+1}^{\Delta t/2}&, ~\textrm{stepsize}~ \frac{\Delta t}{2}, ~ 2 \textrm{ step}
        \end{align*}
        be methods on $[t_i,t_{i+1})$. Then, according to theorem II.23 it holds
        \[
        \begin{aligned}
        y(t_{i+1}) & = y_{i+1}^{\Delta t} + e_p(t_{i+1}) \Delta t^p + O(\Delta t^{p+1}) \\
        y(t_{i+1}) & = y_{i+1}^{\Delta t/2} + e_p(t_{i+1}) \l(\frac{\Delta t}{2}\r)^p + O(\Delta t^{p+1}).
        \end{aligned} \labtag{II.23}
        \]
        Thus, subtracting $(\tfrac{1}{2})^p$ times the first line from the second line, we get
        \[
        y(t_\ipo) - \l(\frac{1}{2}\r)^p y(t_\ipo) = y_\ipo^{\Delta t/2} - \l(\frac{1}{2}\r)^p y_\ipo^{\Delta t} + O(\Delta t^{p+1})
        \]
        and hence
        \[
        y(t_{i+1}) = \frac{y_\ipo^{\Delta t/2} - \l(\frac{1}{2}\r)^p y_\ipo^{\Delta t}}{1-\l(\frac{1}{2}\r)^p} + O(\Delta t^{p+1}).
        \]
        Thus,
        \[
        \tilde y_{i+1}^{\Delta t} = \frac{y_\ipo^{\Delta t/2} - \l(\frac{1}{2}\r)^p y_\ipo^{\Delta t}}{1-\l(\frac{1}{2}\r)^p}
        \]
        is an $O(\Delta t^{p+1})$ approximation of $y(t_{i+1})$.
    \item Using \eqref{II.23} and substracting the second line from the first line, we get
    \[
    O(\Delta t^{p+1}) + \frac{y_\ipo^{\Delta t}-y_\ipo^{\Delta t/2}}{1-\l(\frac{1}{2}\r)^p} = e_p(t_i) \Delta t^p.
    \]
    Setting
    \[
    \textrm{EST}_\ipo^{(\textrm{g})} = \frac{ \norm{ y_{i+1}^{\Delta t}-y_\ipo^{\Delta t/2}}}{1-\l(\frac{1}{2}\r)^p}, \labtag{II.24}
    \]
    for some given tolerance $\textrm{TOL}$ of the global error, we can choose a suitable stepsize, by computing $y_\ipo^{\Delta t}$ and $y_\ipo^{\Delta t/2}$ and then
    \begin{itemize}
        \item increase $\Delta t$, if $\textrm{EST}_\ipo^{(\textrm{g})} \leq \textrm{TOL}$ and
        \item decrease $\Delta t$, if $\textrm{EST}_\ipo^{(\textrm{g})} \geq \textrm{TOL}$.
    \end{itemize}
\end{enumerate}

To derive an estimator of the local error, first observe that $\Delta t$ is not constant, i.e. we write $\Delta t_i$ instead of $\Delta t$. Hence, we want to adjust the time step size according to the local behaviour of the solution. This is why we need an estimate.

Consider a explicit one step method $\Phi$, once with stepsize $\Delta t$ (we denote the outcome by $y_\ipo^{\Delta t}$) and once with stepsize $\Delta t/2$ (we denote the outcome by $\bar y_\ipo^{\Delta t}$. Then,
\[
\bar y_\ipo^{\Delta t} = y_i^{\Delta t} + \frac{\Delta t}{2} \Phi\l(t_i,y_i^{\Delta t},\frac{\Delta t}{2}\r) + \frac{\Delta t}{2} \Phi \l( t_i+\frac{\Delta t}{2}, y_i^{\Delta t}+\Delta t \Phi \l(t_i, y_i^{\Delta t}, \frac{\Delta t}{2}\r), \frac{\Delta t}{2} \r).
\]
By \eqref{II.18}, we know that the local error $\eta(t,\Delta t)=O(\Delta t^{p+1})$ can be written as
\[
\eta(t,\Delta t) = g_{p+1}(t)\Delta t^{p+1} + O(\Delta t^{p+2}).
\]
Thus, it holds
\[
\bar \eta(t_i,\Delta t) = g_{p+1}(t_i) \cdot \l(\frac{1}{2}\r)^{p+1} + g_{p+1} \l(t_i+\frac{\Delta t}{2}\r) \l(\frac{\Delta t}{2}\r)^{p+1} + O(\Delta t^{p+2})
\]
and using a Taylor-expansion in the middle part gives
\[
g_{p+1} \l(t_i+\frac{\Delta t}{2}\r) \l(\frac{\Delta t}{2}\r)^{p+1} = g_{p+1} (t_i) \l(\frac{\Delta t}{2}\r)^{p+1} + O(\Delta t^{p+2}).
\]
Hence, we have
\[
\bar \eta(t_i,\Delta t) = 2 g_{p+1}(t_i) \l(\frac{\Delta t}{2}\r)^{p+1} + O(\Delta t^{p+1}) = \l(\einhalb\r)^p \underbrace{g_{p+1}(t_i) \Delta t^{p+1}}_{=\eta(t,\Delta t)} + O(\Delta t^{p+2}). \tag{$\ast$}
\]
So it holds
\begin{align*}
y(t_\ipo) & = \bar y_\ipo + \bar \eta(t_i,\Delta t) + O(\Delta t^{p+2}) \\
y(t_\ipo) & = y_{i+1} + \eta (t_i,\Delta t) + O(\Delta t^{p+2})
\end{align*}
and plugging in ($\ast$) gives
\begin{align*}
y(t_\ipo) & = \bar y_\ipo + \l(\tfrac{1}{2}\r)^p g_{p+1}(t_i) + \Delta t^{p+1}+O(\Delta t^{p+2})\\
y(t_\ipo) & = y_\ipo + g_{p+1}(t_i) + \Delta t^{p+1}+O(\Delta t^{p+2}).
\end{align*}
We can now define
\[
\textrm{EST}^{(\textrm{l})}_\ipo = \frac{ \norm{\bar y_{i+1}^{\Delta t} - y_\ipo^{\Delta t}}}{\l(\tfrac{1}{2}\r)^p -1}
\]
as an estimate.

This gives the \textbf{Embedding-Stategy for step size control}. Let us assume that we have two one step methods with incremental functions $\Phi_p$ and $\Phi_q$ such that $\Phi_p$ has order of consistency $p$, i.e. $\eta_p(t,\Delta t) = O(\Delta t^{p+1})$ and $\Phi_q$ has order of consistency $q$, i.e. $\eta_q(t,\Delta t) = O(\Delta t^{q+1})$. Without loss of generality, we can assume $q>p$. Then
\begin{align*}
\eta_p(t,\Delta t)
& = y(t) + \Delta t \Phi_p - y(t+\Delta t) + \Delta t \Phi_q - \Delta t\Phi_q \\
& = y(t) + \Delta t \Phi_q - y(t+\Delta t) + \Delta t \l(\Phi_p - \Phi_q\r) \\
& = \eta_q(t,\Delta t) + \Delta t \l(\Phi_p-\Phi_q\r).
\end{align*}
If we now \textbf{assume}\footnote{This is not always given!} that the error $\eta_q$ is negligible with respect to $\eta_p$, we have
\[
|\Phi_p - \Phi_q| \approx O(\Delta t^p)
\]
and so
\[
\eta_p(t,\Delta t)\approx \Delta t(\Phi_p-\Phi_q).
\]
We define
\[
\varepsilon(\Delta t) = \frac{\eta_p(t,\Delta t)}{\Delta t}.
\]
With this, the task of finding a new stepsize $\Delta t_{\textrm{new}}$ is realized by
\[
\varepsilon(\Delta t_{\textrm{new}}) = \varepsilon_{\textrm{goal}}.
\]
With the assumption
\[
\eta_p(t,\Delta t) = c \cdot \Delta t^{p+1}
\]
we obtain
\[
\frac{\varepsilon_{\textrm{goal}}}{\Delta t^p_{\textrm{new}}} \approx \frac{\eta_p(t,\Delta t_{\textrm{new}})}{\Delta t_{\textrm{new}}^{p+1}} \approx c \approx \frac{\eta_p(t,\Delta t)}{\Delta t^{p+1}} \approx \frac{\varepsilon(\Delta t)}{\Delta t^p}
\]
so we can use
\[
\Delta t_{\textrm{new}} = \Delta t \sqrt[p]{\frac{\varepsilon_{\textrm{goal}}}{\varepsilon(\Delta t)}}
\]
as an iteration. Hence, we choose $\Delta t$ courser, if $\varepsilon_{\textrm{goal}}>\varepsilon(\Delta t)$ and else finer.\footnote{Warning! We had $ \eta_p(t,\Delta t) = O(\Delta t^{p+1}), ~ \textrm{for} ~ \Delta t\to0$, so $\eta_p(t,\Delta t)\leq c \Delta t^{p+1}$ only is a reasonable assumption for small $\Delta t$.}

The embedding strategy can be realized in a Runge-Kutta scheme. Consider the RKM given by the Butcher table
\[
\begin{matrix}[c|cccc]
0 & & & &\\
1/2 & 1/2 & & &\\
3/4 & 0 & 3/4  & &\\
1 & 2/9 & 1/3 & 4/9 &\\ \hline
\blau p=2 & \blau 7/24 & \blau 1/4 & \blau 1/3 & \blau 1/8 \\
p=3 & 2/9 & 1/3 & 4/9 &
\end{matrix}
\]
This means
\begin{align*}
r_1 & = f(t_i,y_i) \\
r_2 & = f \l(t_i + \tfrac{\Delta t}{2},~ y_i + \tfrac{\Delta t}{2} r_1\r) \\
r_3 & = f \l(t_i + \tfrac{3}{2}\Delta t,~ y_i + \tfrac{3}{2}\Delta t r_2\r) \\
r_4 & = f \l(t_i + \Delta t,~ y_i + \tfrac{2}{9}\Delta t r_1 + \tfrac{1}{3} \Delta t r_2 + \tfrac{4}{9} \Delta t r_3\r), \\
\end{align*}
and
\begin{align*}
\blau y_{i+1} & \blau =  y_i + \tfrac{7}{24} r_1 + \tfrac{1}{4} r_2  + \tfrac{1}{3}r_3 + \tfrac{1}{8} r_4 \\
y_{i+1} & = y_i + \tfrac{2}{9} r_1 + \tfrac{1}{3} r_2 + \tfrac{4}{9} r_3.
\end{align*}
Hence, it is possible to gain a higher order by simply changing the weights $b_i$. This RKM is due to \textit{Bobacki-Shampine} and can be found in \texttt{MATLAB} as \texttt{ode23}. Note, that the RKM of higher order does not use $r_4$. The Fehlberg-Trick is known by
\[
r_4^i = f\l(t_i+\Delta t, y_i + \Delta t \l( \tfrac{2}{9} r_1^i + \tfrac{1}{3}r_2^i + \tfrac{4}{9} r_3^i\r)\r)=: r_1^{i+1}.
\]
This allows for higher efficiency and is often used in programming.

\subsection{Multistep Methods}
In one step methods, computing $y_{i+1}$ only depends on $y_i$. However, in multistep methods $y_{i+1}$ can also depend on the previous steps $y_{i-1},y_{i-2}$ and so on.

\begin{example}
Let $y_i,y_{i+1}$ be given. Applying the midpoint rule gives
\[
y(t_{i+2}) - y(t_i) = \int_{t_i }^{t_{i+1}} y'(\tau)~ d\tau  = 2 \Delta t f(t_{i+1},y_{i+1})+O(\Delta t^2)
\]
and we thus obtain
\[
y_{i+2} = y_i + 2 \Delta t + f(t_{i+1},y_{i+1}).
\]
\end{example}

\subsubsection*{"Adam's Family"}
Analogue to the last example, we have
\[
y(t_{i+m}) - y(t_{i+m-r}) = \int_{t_{i+m-r}}^{t_{i+m}}f(\tau,y(\tau))~d\tau.
\]
Replacing $f$ by an interpolation polynomial $q$ we obtain
\[
y(t_{i+m}) - y(t_{i+m-r}) = \int_{t_{i+m-r}}^{t_{i+m}}q(\tau)~d\tau.
\]

We now give examples of explicit and implicit multistep methods.

\begin{example}[Adams-Bashorth Method]
Let $r=1$ and $q\in\P_{m-1}$ with
\[
q(t_{i+j}) = f(t_{i+j},y_{i+j})
\]
for $j=0,\dots,m-1$.

For different $m\in\N$, we obtain the methods
\vspace{-1ex}\begin{enumerate}
    \item[$m=1$:] $y_{i+1} = y_i + \Delta t f(t_i,y_i)$ and
    \item[$m=2$:] $y_{i+2} = y_{i+1} + \frac{\Delta t}{2} (3f_{i+1} + f_i)$.
\end{enumerate}
\end{example}

\begin{example}[Adams-Moulton Method]
Let $r=1$ and $q\in\P_m$ with
\[
q(t_{i+j}) = f(t_{i+j},y_{i+j})
\]
for $j=0,\dots,m$.
We obtain the methods
\vspace{-1ex}\begin{enumerate}
    \item[$m=1$:] $y_{i+1} = y_i + \frac{\Delta t}{2} (f_{i+1}+f_i)$,
    \item[$m=2$:] $y_{i+2} = y_{i+1} + \frac{\Delta t}{12} (5 f_{i+2} + 8f_{i+1} + f_i)$ and
    \item[$m=3$:] $y_{i+3} = y_{i+2} + \frac{\Delta t}{24} (9 f_{i+3} + 19 f_{i+2} - 5 f_{i+1} + f_i)$.
\end{enumerate}
\end{example}

\begin{definition}[Multistep Method]
A method for solving the IVP \eqref{II.1} and \eqref{II.2} of the form
\[
\sum_{j=0}^{m }\alpha_j y_{i+j} = \Delta t \Phi(t_i,y_i,\dots,y_{i+m},\Delta t)
\]
with coefficients $\alpha_j\in\R, ~j=0,\dots,m, ~\alpha_m\neq 0$ and given starting values $y_0,\dots,y_{m-1}$ at time steps $t_0,\dots,t_{m-1}$ and incremental function $\Phi: [a,b] \times \R^m \times \R^+ \to \R$ is a \textbf{$m$-step/multistep method}. We call it \textbf{explicit}, if $\Phi$ does not depend on $y_{i+m}$ and \textbf{implicit} otherwise.
\end{definition}

\begin{remark}
\hfill\begin{itemize}
    \item In the following, we confine ourself to linear multistep methods, i.e. to the case
    \[
    \Phi(t_i,y_i,\dots,y_{i+m},\Delta t) = \sum_{j=0}^{m }\beta_j f(t_{i+j},y_{i+j})
    \]
    and claim $|\alpha_0|+|\beta_0|>0$ such that the data at $t_i$ are used for computing $y_{i+m}$.
    \item We need more than one starting values, however often only $y_0$ is given. We thus have to compute the $m-1$ starting values. This is called \textit{starting} of \textit{initialization phase}.
\end{itemize}
\end{remark}

\begin{definition}
A multistep method is called \textbf{consistent of order $p\in\N$} according to the ODE II.1 if for a solution $y$ the \textbf{local discretization error}
\[
\eta(t,\Delta t) = \sum_{j=0}^{m }\alpha_j y(t+j \Delta t) ~ - \Delta t \Phi\Big(t,y(t),y(t+\Delta t),\dots,y(t+m \Delta t), \Delta t\Big)
\]
for $t \in [a,b]$ and $0 < \Delta t < \frac{b-a}{m}$ meets
\[
\eta(t,\Delta t) = O(\Delta t^{p+1}), ~ \Delta t \to 1.
\]
In the case of $p=1$ we call it \textbf{consistent}.
\end{definition}

\begin{theorem}
A linear multistep method has order of consistency $p\in\N$ if and only if it holds
\[
\sum_{j=0}^{m }\alpha_j = 0 \und \sum_{j=0}^{m }a_j j^q = q \cdot \sum_{j=0}^{m} \beta_j j^{q-1} \labtag{II.25}
\]
for $q=1,\dots, p$.
\end{theorem}
\begin{proof}
The expansions
\[
\blau y(t+ j \Delta t) = \sum_{q=0}^{p} \frac{(j \Delta t)^q}{q!} y^{(q)}(t)~+O(\Delta t^{p+1})
\]
and
\[
\gr y'(t+j\Delta t) = \sum_{q=1}^{p} \frac{(j \Delta t)^{q-1}}{(q-1)!} y^{(q)}(t)~ + O(\Delta t^{p+1})
\]
gives
\begin{align*}
\eta(t,\Delta t) & = \sum_{j=0}^{m }\Big(\alpha_j y(t+j \Delta t) - \Delta t \beta_j f\big( t+j \Delta t, y(t+ j \Delta t)\big)\Big) \\
& = \sum_{j=0}^{m }\Big(\alpha_j y(t+j \Delta t) - \Delta t \beta_j y'(t+j \Delta t)\Big) \\
& = \sum_{j=0}^{m }\Bigg(\alpha_j {\blau \bigg(\sum_{q=0}^{p} \frac{(j \Delta t)^q}{q!} y^{(q)}(t)\bigg)} - \Delta t \beta_j{\gr \bigg(\sum_{q=1}^{p} \frac{(j \Delta t)^q}{(q-1)!} y^{(q)}(t)\bigg)} \Bigg) ~+ O(\Delta t^{p+1}) \\
& = {\rot\dots}
\end{align*}
This ends the proof.
\end{proof}

For the further analysis of linear multistep methods, we need the first and second characteristic polynomials. They are given by
\[
\varrho(\xi) = \sum_{j=0 }^{m }\alpha_j \xi^j
\]
and
\[
\sigma(\xi) = \sum_{j=0}^{m} \beta_j \xi^j.
\]
They are also called \textbf{generating polynmials}. We can thus rewrite the conditions \eqref{II.25} in theorem II.26 as
\[
\varrho(1) = \sum_{j=1}^{m} \alpha_j, ~~~ \sigma(1) = \sum_{j=1}^{m} \beta_j, ~~~ \varrho'(\xi) = \sum_{j=1}^{m} j \alpha_j \xi^{j-1}.
\]
For $q=p=1$ we hence have to check
\[
\varrho(1) = 0 \und \varrho'(1) = \sigma(1) \labtag{II.26}
\]
to obtain consistency.

\begin{definition}
A multistep method with starting values $y_j = y(t_j) + O(\Delta t^p), ~\Delta t\to 0$ for $j=0,\dots,m-1$ is convergent of order $p\in \N$ with respect to \eqref{II.1} and \eqref{II.2} if for a time step $\Delta t$ the produced approximation $y_i$ of $y(t_i)$ for $t_i = a+i \Delta t, ~ t_i \in [a,b]$ the \textbf{global error}
\[
e(t_i, \Delta t) = y(t_i) - y_i
\]
for all $t_i$ meets the condition
\[
e(t_i, \Delta t) = O(\Delta t^p), ~ \Delta t\to 0.
\]
If the above conditions also hold for $o(1)$ instead of $O(\Delta t^p)$, we call the method \textbf{convergent}.
\end{definition}

In case of $\beta_j=0$ for all $j$, we have $\Phi = 0$ and thus, the linear multistep method becomes a linear \textbf{homogeneous difference equation}, i.e.
\[
\sum_{j=0}^{m }\alpha_j y_{i+j} = 0. \labtag{II.27}
\]

\begin{theorem}
If the first characteristic polynomial
\[
\varrho(\xi) = \sum_{j=0}^{m} \alpha_j \xi^i \labtag{II.28}
\]
of \eqref{II.27} has only pairwise different roots $\xi_1,\dots,x_m\in\C$, then the solution sequence $\folge{y}$ of \eqref{II.27} is of the form
\[
y_n = \sum_{k=1 }^{m } \gamma_k \xi_k^n, \labtag{II.29}
\]
where $\gamma_k\in\C$ for all $k$.
\end{theorem}
\begin{proof}
We only give a sketch. We need to show
\begin{enumerate}[(a)]
    \item \eqref{II.29} is a solution of \eqref{II.27} and
    \item the solution space has dimension $m$ and $(\xi^n_k)_{n\in\N}$ is a basis.
\end{enumerate}
\begin{enumerate}[{ad }(a).]
    \item By plugging in \eqref{II.29} we obtain
    \[
    \sum_{j=0}^{m }\alpha_j y_{i+j} = \sum_{j=0}^{m} \alpha_j \sum_{k=1}^{m}\gamma_k \xi_k^{i+j} = \sum_{k=1}^{m }\gamma_k\xi_k^i \cdot \underbrace{\sum_{j=0}^{m} \alpha_j \xi_k^j}_{=\varrho(\xi_k)=0} = 0.
    \]
    \item The main idea is to consider that each solution sequence $\folge{y}$ needs starting values $y_0,\dots,y_{m-1}\in\C$. This means the set of starting vectors $\l\{s^{(1)},\dots,s^{(m)}\r\}$ is a basis of $\C^m$. Each starting vector $s=\begin{bmatrix} y_0 & \dots & y_{m-1} \end{bmatrix}^T$ can be written as
    \[
    s=\sum_{j=1}^{m }\gamma_j s^{(j)}.
    \]
    Then, one can show that out of $\folge{s}$ it is possible to build $\folge{y}$.\footnote{This is the long part of the proof.}
\end{enumerate}
\end{proof}

\begin{example}
We consider the ODE
\[
y'(t) = 0,~ y(0) = 0.1 \labtag{II.30}
\]
and apply
\[
y_{i+2} + 4 y_{i+1} - 5 y_i = \Delta t \big( 4f(t_\ipo,y_\ipo) + 2 f(t_i,y_i) \big).
\]
This method is consistent of order $3$ with respect to $y'(t) = f(t,y(t))$. Applying it to \eqref{II.30} gives
\[
y_{i+2} + 4 y_{i+1} - 5y_i = 0.
\]
According to theorem II.28 we have
\[
y_n = \gamma_1 \xi_1^n + \gamma_2 \xi_2^n, \labtag{II.31}
\]
where we have determined $\gamma_1,\gamma_2$ by the initial conditions and $\xi_1$ and $\xi_2$ are the roots of the first characteristic polynomial. We have
\[
\varrho(\xi) = \xi^2 + 4\xi - 5 = (\xi-1)(\xi+5),
\]
so it holds $\xi_1 = 1$ and $\xi_2 = -5$. Setting
\[
y_0 = 0.1 \und y_1=0.1 + \varepsilon,
\]
where $\varepsilon>0$ is a small perturbation, with \eqref{II.31} we obtain
\[
0.1 = \gamma_1 + \gamma_2 \und 0.1 + \varepsilon = \gamma_1 - 5 \gamma_2,
\]
so it holds
\[
\gamma_1 = 0.1 + \frac{\varepsilon}{6} \und \gamma_2 = - \frac{\varepsilon}{6}.
\]
Hence
\[
y_n = 0.1 + \frac{\varepsilon}{6} - \frac{\varepsilon}{6}(-5)^n.
\]
The solution of \eqref{II.30} is $y(t)=0.1$. For a fixed $T$ and $\Delta t=\tfrac{T}{n}$, where $n\in\N$, we have
\[
\limn |y(T)-y_n| = \limn |\frac{\varepsilon}{6} - \frac{\varepsilon}{6} (-5)^n| = \infty.
\]
Even for $y_0=y_1=0.1$ we end up with $|y_{50}|\approx 6.5\cdot 10^{16}$.
\end{example}

\begin{definition}
A multistep method
\[
\sum_{j=0}^{ m }\alpha_j y_{i+j} = \Delta t \Phi(t,y_i,\dots,y_{i+m},\Delta t)
\]
is called \textbf{zero stable}, if the corresponding first characteristic polynomial
\[
\varrho(\xi) = \sum_{j=0}^{ m }\alpha_j \xi^j
\]
meets the \textit{Dahlquist root condition}, i.e. {all roots of the polynomials lie inside of the closed complex unit circle and the roots on the boundary are simple.}
\end{definition}

\begin{remark}
In the beginning of theorem II.26 we had $\varrho(1)=0$. This is why we want the roots on the boundary when they are simple.
\end{remark}

\begin{theorem}
A convergent linear multistep method is necessarily consistent and zero stable.
\end{theorem}
\begin{proof}
We will show
\begin{enumerate}[(i)]
    \item zero stability (by contradiction) and
    \item consistency $\varrho(1) = 0$ and $\varrho'(1)=\sigma(1)$.
\end{enumerate}
\begin{enumerate}[{ad }(i):]
    \item Let us assume, that the method is convergent but not zero stable. Apply the linear multistep method on
    \[
    y'(t) = 0, ~y(0) = 0, ~t\in[0,1]
    \]
    with solution $y(t)=0$. For $i=0,1,\dots$ we have
    \[
    \alpha_m y_{i+m} + \dots + \alpha_0 y_i = 0. \labtag{II.32}
    \]
    Since the method is \textbf{not} zero stable, we have either a simple root $\xi_1$ with $|\xi_1|>1$ or a multiple root $\xi_2$ with $|\xi_2|=1$. Our goal is to generate initial values (meeting the definition for convergence) and show that \eqref{II.32} diverges.

    For the case of $|\xi_1|>1$, we know that according to theorem II.28 a solution of \eqref{II.32} is given as
    \[
    y_n = \sqrt{\Delta t}~\xi_1^n.
    \]
    The first $m$ elements meet
    \[
    \lim_{\Delta t\to0} |y_j - y(t_j)| = \lim_{\Delta t\to0} \l|\sqrt{\Delta t}~\xi_1^j\r| = 0
    \]
    for $j=0,\dots,m-1$ and thus the definition for the convergence. For some fixed $T = n \Delta t \in (0,1)$ with $n>m$ we obtain
    \[
    \lim_{\Delta t\to0} |y_n-y(T)| = \lim_{\Delta t\to0} \l|\sqrt{\Delta t}~\xi^{T/\Delta t}_1\r| = \infty,
    \]
    which is a contradiction to the assumed convergence.

    On the other hand, if $|\xi_2|=1$, since $\xi_2$ is a multiple root it holds
    \[
    0 = \varrho'(\xi_2) = \alpha_1 + 2 \alpha_2 \xi_2 + \dots + m \alpha_m \xi_2^{m-1}.
    \]
    We set
    \[
    y_n = n \cdot \sqrt{\Delta t}~ \xi_2^{n-1},
    \]
    such that
    \[
    \sum_{j=0}^{m} \alpha_j y_{i+j} = \sum_{j=0}^{m }\alpha_j(i+j)\sqrt{\Delta t}~\xi_2^{i+j-1} = \sqrt{\Delta t} ~ i \xi_2^j \sum_{j=0}^{m } j\alpha_j \xi_2^{j-1} = 0
    \]
    is a solution of the difference equation \eqref{II.32}, where the starting values $y_0,\dots,y_{m-1}$ with
    \[
    \lim_{\Delta t\to0} |y_j-y(t_j)| = \lim_{\Delta t\to0} \l| j \cdot \sqrt{\Delta t}\r| = 0
    \]
    meet the definition for convergence. For a fixed $T=n \Delta t\in (0,1)$ with $n\geq m$ we have
    \[
    \lim_{\Delta t\to0} |y_n - y(T)| = \lim_{\Delta t\to0} |n\sqrt{\Delta t} \cdot 1|=\lim_{\Delta t\to0} \l| \frac{T}{\sqrt{\Delta t}}\r| = \infty.
    \]
    This again is a contradiction. Thus, zero stability is necessary for convergence.

    \item For $\varrho(1)=0$ consider the ODE
    \[
    y'(t) =0,~ y(0)=1,~ t\in[0,1]
    \]
    with exact solution $y = 1$. The corresponding difference equation is identical to the one before. For the starting values $y_0,\dots,y_{m-1}$ and for $t = m \cdot \Delta t$ we have
    \begin{align*}
    0 & = \lim_{\Delta t\to 0} \alpha_m (y(m \Delta t)-y_n) \\
    & = \lim_{\Delta t\to 0} \alpha_m(1-y_m) \\
    & = \sum_{j=0}^{m } \alpha_j (1 - y_j) \\
    & = \sum_{j=0}^{m } \alpha_j - \underbrace{\sum_{j=0 }^{m }\alpha_j y_j}_{=0} = \varrho(1).
    \end{align*}

    For showing $\varrho'(1) = \sigma(1)$ we consider
    \[
    y'(t) = 1, ~ y(0) = 0, ~ t\in[0,1].
    \]
    The exact solution is given by $y(t) = t$. This gives the inhomogeneous difference equation
    \[
    \sum_{j=0}^{m }\alpha_j y_{i+j} = \Delta t \sum_{j=0}^{m }\beta_j
    \]
    for all $i$. Since $\varrho(1)=0$, we know that $\xi=1$ is a simple root.\footnote{This follows by the zero stability shown before.} Hence, $\varrho'(1)\neq 0$ and we can set
    \[
    M:= \frac{\sigma(1)}{\varrho'(1)}
    \]
    and
    \[
    y_n:=n\cdot \Delta t \cdot M
    \]
    for all $n\in\N$. The first $m$ elements meet
    \[
    \lim_{\Delta t\to 0 } |y_j - y(j \Delta t)| = \lim_{\Delta t\to 0} |j\cdot \Delta t \cdot (m-1) | = 0
    \]
    for $j=0,\dots,m-1$.\footnote{This means that the starting conditions are met by the first $m$ elements.} Further, it holds
    \begin{align*}
    \sum_{j=0}^{m }\alpha_j y_{i+j} - \Delta t \cdot \sum_{j=0}^{m }\beta_j
    & = \sum_{j=0}^{m }\alpha_j (i+j) \Delta t \underbrace{M}_{=\frac{\sigma(1)}{\varrho'(1)}} - \Delta t \cdot \underbrace{\sum_{j=0}^{m }\beta_j}_{=\sigma(1)} \\
    & = \Delta t \cdot \left( \frac{\sigma(1)}{\varrho(1)} \left( i \cdot \sum_{j=0}^{m }\alpha_j + \sum_{j=0}^{m} j\alpha_j \right) - \sigma(1) \right) \\
    & = 0
    \end{align*}
    which means that $y_n = n \cdot \Delta t \cdot M$ is a solution of the difference equation. Since we have convergence for $T=n\cdot \Delta t \in (0,1)$ fixed, it holds
    \[
    0 = \lim_{\Delta t\to0} |y_n - y(T)| = \lim_{\Delta t \to 0} |n\cdot \Delta t \cdot M - T| = \lim_{\Delta t \to 0} | T - MT| = |T(1-M)|.
    \]
    Hence, since $T$ is nonzero, $1-M$ is zero, so $\sigma(1) = \varrho'(1)$. This ends the proof.
\end{enumerate}
\end{proof}

We now aim to show that
\[
\textrm{zero stability~} + \textrm{~consistency~} \Rightarrow \textrm{~convergence}. \tag{$\ast\ast$}
\]
This leads us to the next definition.

\begin{definition}
A multistep method with incremental function $\Phi$ is \textbf{Lipschitz-continuous in $(t,y(t)$ with Lipschitz constant $L>0$}, if there exists a neighbourhood $U$ of $(t,y(t))$ and a constant $H>0$ such that
\[
\l|\Phi(t_i,u_0,\dots,u_m,\Delta t) - \Phi(t_i, v_0, \dots, v_m, \Delta t)\r| \leq L \cdot \sum_{k=0}^{m } |u_k-v_k|
\]
for all time step sizes $0<\Delta t \leq H$ and all $(t,u_k),(t,v_k)\in U(t,y(t)), ~ k=0,\dots,m$.
\end{definition}

For the proof of our goal we need two technical lemma.

\begin{lemma}[technical]
For $A\in \C^{n,n}$ the sequence of matrices $(A^k)_{k\in\N}$ is bounded if and only if the spectral radius $\varrho(A)\leq1$ and for each eigenvalue $\lambda \in \C$ of $A$ with $|\lambda|=1$ the algebraic and geometric multiplicity are the same.
\end{lemma}
\begin{proof}
We use the Jordan-decomposition $T\inv A T = J = \textrm{diag}(J_1,\dots,J_r)$.\footnote{It states that for each matrix $A\in\C^{n,n}$ there exists a nonsingular transformation $T\in\C^{n,n}$ such that $T\inv A T = J$, where $J= \textrm{diag}(J_1,\dots,J_r)$, where $J_i$ a is a Jordan block.} In each $J_j$ with $m_j\geq 2$, we know that $|\lambda_j|<1$. We now find a $\varepsilon>0$ such that
\[
|\lambda_j|+\varepsilon < 1.
\]
an then we multiply $J$ by $D=\textrm{diag}(1,\varepsilon,\varepsilon^2, \dots, \varepsilon^{k-1})$ and obtain
\[
\tilde J = D\inv J D.
\]
We arrive at the Jordan blocks
\[
\tilde J_j = \begin{bmatrix}
\lambda_j & \varepsilon & & 0 \\
& \ddots & \ddots & \\
& & \lambda_j & \varepsilon \\
0 & & & \lambda_j
\end{bmatrix}.
\]
For $k\in\N$ it holds
\[
\norm{\tilde J^k} \leq \norm{\tilde J}^k \leq \varrho(\tilde J)^k \leq 1. \tag{$\ast$}
\]
Setting $S= D\inv T$, we have $A^k = (S \tilde J S\inv)^k = (S\tilde J^k S\inv)$. Using $(\ast)$ we have
\[
\norm{A^k}_{\infty} \leq \norm{S}_\infty \norm{S\inv}_\infty \norm{\tilde J^k}_{\infty} \leq \kappa_\infty(S) < \infty
\]
for all $k\in \N$.

The reverse direction is much shorter and is an exercise.
\end{proof}

\begin{lemma}[discrete Gronwall-Lemma]
Let $\Delta t_0,\dots , \Delta t_{r-1}\in \R_{>0}$ and $\delta,\gamma\in\R_{\geq 0}$ be given. If the numbers $e_0,\dots,e_r\in\R$ meet
\[
|e_0|\leq \delta \und |e_\ell| \leq \delta + \gamma \sum_{j=1}^{\ell-1} \Delta t_j ~|e_j|
\]
for $\ell=1,\dots,r$, then it holds
\[
|e_\ell| \leq \delta \cdot \exp \l(\gamma \sum_{j=0}^{\ell-1}\Delta t_{j}\r)
\]
for $l=0,\dots,r$.
\end{lemma}
\begin{proof}
No proof, since the proof is similar to the known Gronwall-Lemma.
\end{proof}

\begin{theorem}
Let the multistep method
\[
\sum_{j=0}^{m }\alpha_j y_{i+j} = \Delta t \Phi(t_i, y_i, \dots, y_{i+m},\Delta t) \labtag{II.33}
\]
be Lipschitz and zero stable. Then there exists a stepsize bound $H>0$: $0<\Delta t = \frac{b-a}{n} \leq H$ such that
\[
\max_{j=0,\dots,n} |e(t_j,\Delta t)| \leq K \l( \max_{k=0,\dots,m-1} |e(t_k,\Delta t)| + \max_{a\leq t \leq b - m \Delta t} \frac{\eta(t,\Delta t)}{\Delta t} \r) \labtag{II.34}
\]
with a constant $K>0$ being independent of the Lipschitz-constant.
\end{theorem}
\begin{proof}
For the left hand side of \eqref{II.34}, let us define
\[
e := \begin{bmatrix} e_0 & \dots & e_n  \end{bmatrix}^T \in \R^{n+1},\footnote{Note that $e$ is a vector.}
\]
where $e_i = e(t_i,\Delta t),~i=0,\dots,m$. Furthermore, for each $j=0,\dots,n-(m-1)$ we define
\[
e_j := \begin{bmatrix} e_j & \dots & e_{j+m-1} \end{bmatrix}^T\in\R^m
\]
such that we can express the total error as
\[
\norm{e}_\infty = {\rot\dots}
\]
Furthermore, we set
\[
\eta_i = \eta(t_i,\Delta t), ~i=0,1,\dots
\]
and obtain
\begin{align*}
\sum_{j=0}^{m }\alpha_j e_{i+j} & = \sum_{j=0}^{m }\alpha_j \l( y(t_{i+j}) - y_{i+j} \r) \\
& = \Delta t \underbrace{\l( \Phi\Big(t_i,y(t_i),\dots,y(t_{i+m}),\Delta t\Big) - \Phi\Big(t_i, y_i, \dots, y_{i+m},\Delta t\Big) \r)}_{=:\mu_i} + \eta_i. \labtag{II.35}
\end{align*}
Without loss of generality, we set $\alpha_m=1$. With \eqref{II.35} we obtain
\[
e_{i+m} = \sum_{j=0}^{m-1} \alpha_j e_j + \mu_i + \eta_i.
\]
Finally, we arrive at
\[
\begin{bmatrix} e_{i+1} \\ \vdots \\ e_{i+m} \end{bmatrix}
= \begin{bmatrix}
0 & 1 & & \\
& \ddots & \ddots& \\
& & 0 & 1 \\
-\alpha_0 & \dots & \dots & -\alpha_{m-1}
\end{bmatrix}  \begin{bmatrix} e_i \\ \vdots \\ \vdots \\ e_{i+m-1} \end{bmatrix} + \begin{bmatrix} 0 \\ \vdots \\ 0 \\ \mu_i+\eta_i \end{bmatrix},
\]
or
\[
e_{i+1} = A e_i + f_i.
\]
We thus can rewrite
\[
e_i = A e_{i-1} f_{i-1} = \dots = A^i e_0 + \sum_{k=0}^{i-1 } A^{i-1-k} f_k. \labtag{II.36}
\]
This means that we can express the error vector $e_i$ in terms of the initialization error. One can show that
\[
p(\lambda)=\det(A-\lambda I) = - \sum_{j=0}^{m-1} \alpha_j \lambda^j = - \varrho(\lambda).
\]
Due to the zero stability we know that for all eigenvalues $\lambda$ it holds $|\lambda|\leq 1$ and the eigenvalues with $|\lambda|=1$ are simple. We can thus apply lemma II.32 such that $(A^k)_{k\in\N}$ is bounded, i.e.
\[
\norm{A^p}_{\infty} \leq C
\]
for all $p\in \N$ with a $C\geq 1$. We thus rewrite \eqref{II.36} as
\[
\norm{e_i}_\infty \leq C \l( \norm{e_0}_\infty  + \sum_{k=0}^{i-1} \norm{f_k}_\infty  \r). \labtag{II.37}
\]
We still need to handle $\norm{f_k}_\infty$. Our idea is to use the Lipschitz condition
\[
\norm{\mu_k} \leq \Delta t \cdot L \cdot \sum_{j=0}^{m} |y(t_{k+j}) - y_{k+j}| = \Delta t \cdot L \cdot \sum_{j=0}^{m }|e_{k+j}|.
\]
Using that, we obtain
\begin{align*}
\norm{f_k}_\infty
& = |\mu_k+\eta_k| \\
& \leq |\mu_k| + |\eta_k| \\
& \leq |\eta_k| + \Delta t\cdot L \cdot \sum_{j=0}^{m }|e_{k+j}|  \\
& \leq n\cdot \max_{j=0,\dots,n-m} |\eta_j| + \Delta t \cdot L \cdot m \cdot \norm{e_k}_\infty + \Delta t \norm{e_{k+1}}_\infty.
\end{align*}
It further holds
\[
\sum_{k=0}^{i-1}\norm{f_k}_\infty \leq n \cdot \max_{j=0,\dots,n-m} |\eta_j| + \Delta t \cdot L(m+1) \cdot \sum_{k=0}^{i-1} \norm{e_k}_\infty + \Delta t \cdot L \cdot \norm{e_i}_\infty.
\]
Inserting these two results into \eqref{II.37} gives
\begin{align*}
& \norm{e_i}_\infty \\
& \leq C \l( \norm{e_0}_\infty + n \cdot \max_{j=0,\dots,n-m}|\eta_j| +  \Delta t\cdot(m+1) \cdot L \sum_{k=0}^{i-1} \norm{e_k}_\infty + \Delta t \cdot L \cdot \norm{e_i}_\infty \r).
\end{align*}
We thus have
\[
(1-C\cdot \Delta t\cdot L )\norm{e_i}_\infty \leq C \l( \norm{e_0}_\infty + n \cdot \max_{j=0,\dots,n-m} |\eta_j| + \Delta t \cdot (m+1)\cdot L \sum_{k=0}^{i-1}\norm{e_k}_\infty \r).
\]
Since we bound the time step $\Delta t$ by a constant $H<\tfrac{1}{CL}$, we have
\[
(1- \Delta t \cdot C \cdot L ) \geq 1 - HCL > 0.
\]
We thus have
\[
\norm{e_i}_\infty \leq \delta + \gamma \cdot \sum_{k=0}^{i-1} \Delta t \norm{e_k}_\infty
\]
for $i=0,\dots,n-m$, where
\[
\delta:= \frac{C}{1-HCL} \l( \norm{e_0}_\infty + n \cdot \max_{j} |\eta_j| \r) \geq 0
\]
and
\[
\gamma:= \frac{(m+1)\cdot L \cdot C}{1-HCL} \geq 0.
\]
With Gronwall's lemma (lemma II.33), we obtain
\[
\norm{e_i}_\infty \leq \delta \cdot \exp \l( \gamma \cdot \sum_{k=0}^{i-1}\Delta t \r) \leq \delta\cdot \exp(\gamma(b-a)).
\]
Setting
\[
\tilde k:= \frac{C}{1-HCL}\cdot \exp(\gamma(b-a)) \textrm{~and~} n = \frac{b-a}{\Delta t},
\]
we arrive at
\begin{align*}
|e_i| \leq \norm{e_i}_\infty
& \leq \tilde k \l( \norm{e_0}_\infty + n\cdot \max_{j=0,\dots,n-m}|\eta_j| \r) \\
& \leq \tilde k \l( \norm{e_0}_\infty + n \cdot \max_{a\leq t \leq b-m\Delta t} |\eta(t,\Delta t)| \r).
\end{align*}
Thus, it holds
\[
|e_i| \leq k \cdot \l( \max_{j=0,\dots,m-1} |y(t_i)-y_i| + n\cdot \max_{a\leq t \leq b-m \Delta t} \frac{n(t,\Delta t)}{\Delta t} \r)
\]
for $i=0,\dots,n-(m-1)$ with
\[
k:= \max \l\{\tilde k, \tilde k \cdot (b-a)\r\} \geq 0.
\]
Because of
\[
\max_{j=n-(m-1),\dots,n}|y(t_j)-y_j| = \max_{n-(m-1),\dots,n}|e_j| \leq \norm{e_{n-(m-1)}}_\infty
\]
we obtain
\[
\max_{i=0,\dots,n} |y(t_i)-y_i| \leq k \l( \max_{k=0,\dots,m-1} |e(t_k,\Delta t)| + \max_{a\leq t \leq b-m \Delta t} \frac{|\eta(t,\Delta t)|}{\Delta t} \r),
\]
what ends the proof.
\end{proof}

\begin{remark}\hfill
\begin{itemize}
    \item It holds $e(t_i,\Delta t)=y(t_i)-y_i$.
    \item Theorem II.34 in short states, that for a Lipschitz-continuous multistep method it holds $(\ast\ast)$.
    \item Starting values need at least the order of consistency of the method.
    \item Compare $\frac{\eta(t_i,\Delta t)}{\Delta t}$ with the consistency definition.
\end{itemize}
\end{remark}

\subsubsection*{BDF Methods (Backwards Differentiation Formula Methods)}

\begin{definition}[BDF-scheme]
For $t_{i+j} = t_i + j \Delta t, ~ j=0,\dots,m$ and nodes $(t_i,y_i), \dots, (t_{i+m-1},y_{i+m-1}) \in \R^2$ let $p\in\P_m$ be given via
\[
p(t_{i+j})=y_{i+j} \und p'(t_{i+m}) = f(t_{i+m},p(t_{i+m})) \labtag{II.38}
\]
for $j=0,\dots,m-1$. Then, the BDF scheme computes the approximation $y_{i+m}$ at $t_{i+m}$ by
\[
y_{i+m} = p(t_{i+m}). \labtag{II.39}
\]
\end{definition}

\begin{example}[BDF(2)-method]
For given two values $y_i,y_{i+1}$ we consider $p\in\P_2$ such that
\[
p(t_i) = y_i, ~ p(t_{i+1}) = y_{i+1} \und p'(t_{i+2}) = f(t_{i+2},p(t_{i+2})). \tag{$\ast$}
\]
Using an expansion, we get the equations
\[
p(t_{i+1}) = p(t_{i+2}) + \Delta t p'(t_{i+2}) + \frac{\Delta t^2}{2}p''(t_{i+2})
\]
and
\[
p(t_i) = p(t_{i+2}) + 2 \Delta t p'(t_{i+2}) + 2 \Delta t^2 p''(t_{i+2}).
\]
Substracting the first equation from the second four times, we obtain
\[
p(t_i) - 4 p(t_{i+1}) = -3 p(t_i) - 2 \Delta t p'(t_{i+2}) + 0.
\]
The conditions ($\ast$) give
\[
\frac{y_i - 4 y_{i+2} + 3 y_{i+2}}{2 \Delta t} = p'(t_{i+2}) = f(t_{i+2},y_{i+2}).
\]
Rearranging leads to
\begin{align*}
\frac{3 y_{i+2} - 4 y_{i+1} + y_i }{2 \Delta t}
& = \frac{3 p(t_{i+2}) - 4 p(t_{i+1}) + p(t_i)}{2 \Delta t} \\
& = p'(t_{i+2}) \\
& = f(t_{i+2},p(t_{i+2})) \approx y'(t_{i+2})
\end{align*}
\end{example}

\begin{remark}\hfill
    \begin{itemize}
        \item Due to \eqref{II.38} and \eqref{II.39} all BDF methods are implicit.
        \item It can be shown, that BDF(2) is zero stable and consistent of order $p=2$ and also convergent of order $p=2$. The proof relies on Taylor expansion.
    \end{itemize}
\end{remark}

\subsubsection*{Stiff ODEs}
\begin{example}[mass-spring system]
Let $m=25 \textrm{kg}$, spring sonstant $D=5000 \textrm{N/m}$, damping $\sigma = 4000 \textrm{kg/s}$ and consider
\[
m u''(t) + \sigma u'(t) + D u(t) = 0, ~ u(0) = 0.05 \textrm{m}, ~ u'(0) = 10 \textrm{m/s}.
\]
The roots of the characteristic polynomial are given by
\[
\lambda_{1,2} = -80 \pm \sqrt{6200},
\]
i.e.
\[
\lambda_1 \approx -158.74 \und \lambda_2 \approx -1.26.
\]
The solution is approximated by
\[
u(t) \approx - 0.064 \cdot \exp(-158.74 t) + 0.114 \exp(-1.26 t). \footnote{We ignore the units.}
\]
Remember that we wanted $|R(\Delta t\cdot\lambda)|<1$. But since the first summand is very small, the second summand is dominant and the first summand is almos negligible. Regarding $\lambda_2$, the stability is easy to "buy" by choosing $\Delta t = 0.5$, but for $\lambda_1$ we need a much smaller $\Delta t$, despite the fact that its influence is much smaller.

The dominance of the $\lambda_2$ in combination with the small stepsize needed for the $\lambda_1$ part is our problem. Due to stability reasons, the whole interval needs fine discretization. Hence, the eigenvalue with the largest absolute value is important for the stability.

In the literature, one can find different measures for the stiffness. A popular one is
\[
S:=\frac{\max_j |\textrm{Re}(\lambda_j)|}{\min_j |\textrm{Re}(\lambda_j)|}.
\]
In our example, we have $S\approx 100$.

Solving our problem in \texttt{MATLAB} needs 1957 time steps when using \texttt{ode45}\footnote{This is a lot.} and 72 time steps when using \texttt{ode23s}. The reason behind this is that \texttt{ode23s} uses a more stable method, the \emph{Rosenbrock method}.
\end{example}

\subsubsection*{Rosenbrock-Method}
Consider
\[
y' = f(y), ~J:=f'(y_n), ~ y'=f(f(y_n)+J(y-y_n)) + (f(y)-f(y_n)) - J(y-y_n),
\]
i.e. a linearization of $f$ at $y_n$ and assume
\[
y' = Jy(t) + f(y(t)) - Jy(t).
\]
We consider the slopes $r_j$\footnote{Compare with the Runge-Kutta-scheme.}
\[
r_{\gr j} = J \l( y + \Delta t \sum_{\nu=1}^{\gr j}\tilde{a}_{j\nu} r_\nu \r) + \l( f\l(y+\Delta t \sum_{\nu=}^{j-1} a_{j\nu} r_\nu\r) - J \l( y+ \Delta t \sum_{\nu=1}^{j-1} a_{j\nu} r_\nu \r) \r).
\]
Hence, $r_{\gr j}$ appears only linear on the right hand side. Further, $r_1,\dots,r_{j-1}$ appear only in the nonlinear part of $f$. This means that we can compute $r_j$ via a linear system:
\begin{enumerate}[(i)]
    \item Compute $J= f(y_n)$.
    \item Compute
    \[
    (I-\Delta t \tilde a_{jj} J ) r_{j} = \Delta t \sum_{\nu=1}^{j-1} \l( \tilde a_{j\nu} - a_{j\nu} \r) J \cdot r_0  ~ + ~ f \l( y_n + \Delta t \sum_{\nu=1}^{j-1} a_{j\nu} r_\nu \r).
    \]
    \item Compute
    \[
    y_{n+1} = y_n + \Delta t \sum_{ j=1}^{ s} b_j r_j.
    \]
\end{enumerate}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage
\section{Numerics of Boundary Value Problems (BVPs)}
In IVPs, we had initial conditions. In BVPs, there are conditions on the solution, usually on the boundary of $I=[x_0, x_{\textrm{end}}]\subseteq \R$.

\includepdf[pages={1-9},nup=2x3]{BVP_slides.pdf}

\begin{example}
Consider the ODE $-u''(x) = u(x)$ with boundary conditions $u(0)=0, u(\pi/2)=1$. Then the unique solution is given by $u(x) = \sin(x)$. When choosing the boundary conditions $u(0) = 0, u(\pi)=0$, we get infinitely many solutions of the form $u(x)=c_1 \cdot \sin(x)$. Furthermore, choosing $u(0)=0, u(\pi)=1$ there is no solution.
\end{example}

In general, for each BVP
\[
y' = f(x,y), ~ R(y(a),y(b)) = 0
\]
we can give an IVP
\[
y' = f(x,y), ~ y(a) = s
\]
with still unknown parameter $s$. If these initial conditions give unique solutions
\[
w(x,s):= y(x;a,s),
\]
then there exists a unique function
\[
F(x) = R(s,w(b,s)).
\]
Thus, the solution is uniquely solvable if $F(s)=0$ has exactly one solution.

Let $(u_i)_{i\in\N}$ be a sequence of functions such that
\[
\norm{u_i - y}_\infty = \delta_i, ~ \delta_i \to 0
\]
and instead of $y(a)=s^*$ we have $u_i(a)=s^* + \alpha_i$ with $\alpha_i \to 0$ and
\[
u_i(b) = w(b,s^*) + \beta_j, ~ \beta_j \to 0.
\]
In other words, $(u_i)_{i\in\N}$ is a sequence of functions approaching the solution. We use the expansion
\[
\frac{F(s^* + \alpha_i)}{\norm{\alpha_i}} = \l( F(s^*) + F_s(s^*) \alpha_i + O(\alpha_i^2)\r) \cdot \einsdurch{\norm{\alpha_i}}.
\]
Since $F(s^*)=0$, we only need to consider the part $F_s(s^*) \cdot \alpha_i$. In other words,
\[
\frac{F(s^* + \alpha_i)}{\norm{\alpha_i}}=0 \Rightarrow F_s(s^*) \alpha_i \to 0.
\]
If $F_s$ is continuous, then $F_s(s^*)\cdot \alpha = 0$. But since
\[
F(s) = R(s,w(b,s))
\]
it holds
\[
F_s(s^*) \cdot \alpha = \partiell{R(s,w(b,s))}{y(a)} \alpha + \partiell{R(s,w(b,s))}{y(b)} \cdot \partiell{w(b,s)}{s}\cdot\alpha.
\]
If the solution is not local unique, then $F_s$ is singular and there exists a direction $\alpha$ for varying the $s^*$ such that
\[
\partiell{R(s,w(b,s))}{y(a)}\cdot \alpha + \partiell{R(s,w(b,s))}{y(b)} \cdot \underbrace{\partiell{w(b,s)}{s} \cdot \alpha}_{=:\beta} = 0.
\]

\subsection{Finite Elements}
We consider
\[
\begin{aligned}
-u''(x) + b(x) u'(x) + c(x) u(x) & = f(x) \\
0< x < 1, ~ u(0) = u(1) = 0
\end{aligned} \labtag{III.4}
\]
and expect that there exists a solution
\[
u \in C^2(0,1) \cap C^0([0,1]).
\]
Such a solution is called \emph{classical solution}.

We seek for a minimization problem, but it can happen that we can find a suitable minimization, but it has solution or it has only a solution by assuming that the solution is suitable smooth.

Let
\[
X:= \l\{ v\in C^1(0,1) \cap C^0([0,1]) ~|~ v(0) = v(1) = 0 \r\}.
\]
Multiplying \eqref{III.4} by some $v\in X$ gives
\[
\int_0^1 \l({\gr - u''} + b  u' + c  u \r)  {\gr v} ~ dx =  \int_{0}^{1}f v ~ dx.
\]
With integration by parts on $\gr -u''  v$ gives
\[
\underbrace{\int_{0}^{1}{\gr - u'' v} + b u' v + c u v}_{=:a(u,v)}~ dx = \underbrace{-u'(x) v(x)\Big|^1_0}_{=0} + \underbrace{\int_{0}^{1}u' v' + b  u'  v + c u  v~dx}_{=a(u,v)} = \underbrace{\int_{0}^{1}f v~dx}_{=:f(v)=\scal{f}{v}}
\]
for all $v\in X$.

\setcounter{theorem}{4}
\begin{definition}
The closure (\emph{Abschluss}) of a subset $A$ of $V$ with respect to $V$ is denoted by $\overline{A}^V$. By
\[
\supp v := \ol{\l\{x\in (0,1)~|~ v(x) \neq 0\r\}}^{\R}
\]
we define the \textbf{support}.
\end{definition}

\begin{remark}\hfill
\begin{itemize}
    \item For an example, let $(X,\norm{\cdot}_X)$ be some normed space and $A\subseteq X$. Then $x \in \ol{A}^{\norm{\cdot}_X}$ if there exists $\folge{x}\subseteq A$ such that
    \[
    \limn \norm{x_n-x}_X = 0.
    \]
    \item We are interested in functions with small support.
\end{itemize}
\end{remark}

\begin{definition}
We define
\[
C_0^\infty(0,1):= \l\{ v \in C^\infty ~|~ \supp v \subseteq (0,1)\r\}
\]
and
\[
L^2(0,1):= \l\{ v: [0,1] \to \R ~|~ \int_{0}^{1}v^2(x)~dx < \infty \r\}
\]
and equip $L^2(0,1)$ with the norm
\[
\norm{u}^2_{L^2(0,1)} = \int_{0}^{1}u^2(x)~dx
\]
and inner product (\emph{Skalarprodukt})
\[
\scal{u}{v}_{L^2(0,1)} = \int_{0}^{1}uv~dx
\]
\end{definition}

\begin{remark}
    \item The sets $L^p(\Omega)$ (equipped with the respective norm) are for $1\leq p < \infty$ are separable Banach spaces.
\end{remark}

\begin{definition}
We call $w\in L^2(0,1)$ the \textbf{weak derivative} of $u \in L^2(0,1)$ if for all $v\in C^\infty_0(0,1)$ it holds
\[
\int_{0}^{1}w v ~ dx = - \int_{0}^{1}u v'~ dx.
\]
In this case, we write $w=u'$.
\end{definition}

\begin{definition}
We call the space
\[
H^1(0,1) := \l\{ v\in L^2(0,1) ~|~ \exists v' \in L^2(0,1) \r\}
\]
\textbf{sobolev space of the functions on $(0,1)$}. Further, we set
\[
H_0^1(0,1) = \ol{C^\infty_0(0,1)}^{H^1(0,1)}.
\]
\end{definition}

But what does $\ol{C_0^\infty(0,1)}^{H^1(0,1)}$ mean?

\begin{remark}\hfill
\begin{enumerate}[(1)]
    \item The set $H_0^1(0,1)$ is the closure of $C_0^\infty(0,1)$ with respect to $H^1(0,1)$.
    \item It holds $u\in H_0^1(0,1)$ if and only if there exists a sequence of functions $\folge{u}\in C^1_0(0,1)$ such that
    \[
    \norm{u-u_n}_{H^1(0,1)} \overset{n\to\infty}{\longrightarrow} 0
    \]
\end{enumerate}
\end{remark}

\begin{definition}
The solution of the problem
\[
\textrm{"Find } u \in H^1_0(0,1) \textrm{ such that } a(u,v) = f(v) \textrm{ holds for all } v\in H^1_0(0,1) \textrm{."} \labtag{III.5}
\]
is called \textbf{weak solution} of (III.1)
\end{definition}

\subsubsection*{Minimization Problem}
Consider
\[
u''(x) = f(x) ~\forall x \in(0,1), ~ u(0)=u(1)=0.
\]
Analogue to \eqref{III.5} we have
\[
a(u,v) = f(v) \und a(u,v) = \int_{0}^{1}u'v'~dx.
\]
This gives the minimization problem
\[
\textrm{"Find } u \in V \textrm{ such that } (u,v) = \scal{f}{v} \textrm{ holds for all } v \in V \textrm{".}
\]
For this setting, we consider
\[
J(v) = \frac{1}{2} (u,v) - \scal{f}{v}.
\]
Let us assume that we have a minimizer $u$, i.e.
\[
J(u) = \min_{v \in V} J(v). \tag{$\ast$}
\]
That means
\[
J(u) \leq J(u+tw) ~\forall w\in V, \forall t\in [0,1].
\]
We have
\begin{align*}
J(u+tw)
& = \einhalb \l( (u+tw,u+tw) \r) - \scal{f}{u+tw} \\
& = \einhalb \l( (u,u) + 2 t (u,w)  + t^2 (w,w) \r) - \scal f u - t \scal f w \\
& = J(u) + \frac{t^2}{2} \norm{w}^2 + t (u,w)  - t \scal f w.
\end{align*}
Thus, ($\ast$) is equivalent to the problem
\[
\textrm{"Find } u \in V \textrm{ such that } \frac{t}{2} \norm{w}^2 + (u,w) - \scal f w \geq 0 \textrm{ for all }  w \in V \textrm{ and } t \in [0,1] \textrm{."}
\]
which is equivalent to finding $u \in V$ such that
\[
(u,w) \geq \scal f w
\]
for all $w \in V$.\footnote{Since $w\in V$ implies $-w\in V$ we have an equality here.}

\subsubsection*{Ritz-Galerkin Problem}
We want to solve the problem
\[
\textrm{"Find } u \in V \textrm{ such that } a(u,v) = f(v) \textrm{ for all } v\in V  \textrm{."} \labtag{III.6}
\]
with a bilinear form $a(\cdot,\cdot): V\times V\to \R$ and a linear functional $f: V \to \R$ in some Hilbert space $V$.

Choose $V_n \subseteq V$ with $\dim(V_n)=n <\infty$ and solve
\[
\textrm{"Find } u_n \in V_n \textrm{ such that } a(u_n,v) =f(v) \textrm{ for all } v \in V_n \textrm{."} \labtag{III.7}
\]
Let $\{\varphi_i\}_{i\in\N}$ be a basis of $V_n$. Then \eqref{III.7} can be rewritten as
\[
\textrm{"Find } u_n \in V_n \textrm{ such that } a(u_n,\varphi_i) = f(\varphi_i) \textrm{ for all } i\in \N \textrm{."}
\]
with
\[
u_n = \sum_{i=1}^{n} \ol{u_i} \varphi_i.
\]
We arrive at a linear system for the coefficients $\ol{u_i}$, i.e.
\[
a(u_n,\varphi_i) = a\l(\sum_{j=1}^{n}\ol{u_j} \varphi_j(x), \varphi_i(x) \r) = \sum_{j=1}^{n} \ol{u_j} a(\varphi_j(x),\varphi_i(x)) = f(\varphi_i)
\]
or
\[
A_n U_n = F_n,
\]
where
\[
U_n := \begin{bmatrix} \ol{u_1} & \dots & \ol{u_n} \end{bmatrix}^T,~~A_n := [a_{ij}], ~ a_{ij}:= a(\varphi_j,\varphi_i) ~ \und F_n := \begin{bmatrix} f(\varphi_1) & \dots & f(\varphi_n) \end{bmatrix}^T.
\]

The next theorem gives the answer to the question whether this linear system has always a solution.

\begin{theorem}
Let $V$ be a normed space, $V_n \subseteq V$ a subspace with $\dim V_n = n < \infty$ and $a(\cdot,\cdot): V\times V \to \R$ an $V$-elliptic (coercive) bilinear form, i.e. there exists a $\gamma>0$ such that
\[
a(v,v) \geq \gamma \norm{v}^2_V
\]
for all $v\in V$ and $f: V \to \R$ linear and continuous, i.e. there exists $k>0$ such that
\[
|f(v)| \leq k \norm{v}_V
\]
for all $v\in V$. Then,
\begin{enumerate}[(i)]
    \item the linear system $A_n U_n = F_n$ has a unique solution and
    \item it holds $\norm{u_n}_V \leq \tfrac{k}{\gamma}$.
\end{enumerate}
\end{theorem}
\begin{proof}
Let $W_n = \begin{bmatrix} \ol{w_1} & \dots & \ol{w_n} \end{bmatrix}^T \neq 0$ and
\[
w = \sum_{i=1}^{n}\ol{w_i}\varphi_i.
\]
Then
\[
A_n W_n W_n = a(w,w) \geq \gamma \norm{w}^2_V > 0,
\]
hence $A_n W_n \neq 0$ for all $W_n \neq 0$. This shows that $A_n$ is invertible and hence a unique solution exists. Further, with this, it holds
\[
\gamma \norm{u_n}^2_V \leq a(u_n,u_n) = f(u_n) \leq k \norm{u_n}_V.
\]
\end{proof}

Now, we are concerned with the relation between the finite dimensional and the variational problem and how it is possible to construct our scheme such that $A_n U_n = F_n$ is numerically "easy" to solve.

\begin{theorem}(Lemma of Cea)
Let $V$ be a Hilbert space, $V_n\subseteq V$ be a subspace with $\dim V_n = n < \infty$ and
\[
a(\cdot,\cdot): V\times V \to \R
\]
be a $V$-elliptic and continuous bilinear form, where continuity means
\[
\exists \Gamma > 0: ~ a(u,v) \leq \Gamma \norm{u}_V \norm{v}_V ~ \forall u,v \in V.
\]
Then, it holds
\[
\norm{u-u_n}_V \leq \frac{\Gamma}{\gamma} \inf_{v\in V_n}\norm{v-u}_V.
\]
\end{theorem}
\begin{proof}
Since $V_n \subseteq V$, we have
\[
u \in V: ~ a(u,v_n) = f(v_n) ~ \forall v_n \in V_n
\]
and
\[
u_n \in V_n: ~ A(u_n,v_n) = f(v_n) ~ \forall v_n \in V_n.
\]
Substracting the second equation from the first, we get
\[
a(u-u_n, v_n) =0 ~ \forall v_n \in V_n.
\]
This is called Galerkin-orthogonality. For all $u_n\in V_n$, we have
\begin{align*}
\gamma \norm{u-u_n}_V^2
& \leq a(u-u_n,u-u_n) \\
& = a(u-u_n,u-v+v-u_n) \\
& = a(u-u_n,u-v) + \underbrace{a(u-u_n,v-u_n)}_{=0} \\
& \leq \Gamma \norm{u-u_n}_V \norm{u-v}_V.
\end{align*}
This completes the proof.
\end{proof}

\begin{remark}
If $V_n$ has polynomials of degree $k$, then the approximation quality of the interpolation of $u$ in $\P_k$ can be transferred tothe approximation quality of $u-u_n$ in $V_n$.
\end{remark}

\subsubsection*{Choice of basis functions}
Consider
\[
-u'' = f ~ \textrm{in}~ (0,1), ~ u(0)=u(1)=0
\]
with
\[
a(u,v)=\int_{0}^{1}u'v'~dx
\]
and
\[
a_{ij} = a(\varphi_i,\varphi_j) = \int_{0}^{1}\varphi_i'\varphi_j'~dx.
\]
For instance, choose basis functions
\[
\varphi_j(x) = \sin(j\pi x), ~ j=1,\dots,n \footnote{One can show that this is a basis.}
\]
and set the meshsize $h=\tfrac{1}{n}$ and
\[
V_n = \spann \{\varphi_i\}_{i=1}^n = \l\{ v ~|~ v(x) = \sum_{j=1}^{n}c_j \varphi_j(x) \r\} \subseteq V
\]
It can then be shown, that
\[
\lim_{h\to 0} \l( \inf_{v\in V_n} \norm{u-v} \r) = 0
\]
for any given $u\in V$. We then have
\[
a(\varphi_i,\varphi_j) = \begin{cases} \pi^2 \cdot \frac{j^2}{2} & \textrm{ if } i=j \\ 0 & \textrm{ if } i\neq j \end{cases}.
\]
For the right hand side, we obtain
\[
F_i = \int_{0}^{1} f(x) \varphi_i(x)~ dx
\]
and it can then be shown, that the solution
\[
\sum_{j=1}^n a(\varphi_i,\varphi_j) \cdot \ol{u_j} = f(\varphi_j) ~ \forall i=1,\dots,n
\]
is
\[
\ol{u_j} = \frac{2 F_i}{\pi^2 j^2}.
\]
The Galerkin approximation is then
\[
u_n(x) = \sum_{j=1}^{n} \ol{u_j} \sin(j\pi x).
\]

But having an explicit solution is not given in general. This case is exceptional, since the basis functions $u_j$ are the eigenfunctions of the differential operator.

Consider the same problem as before, but now choose
\[
V_n := \spann \l\{ \einsdurch{i} x^i \r\}_{i=1}^n.
\]
The Galerkin approximation is then
\[
u_n(x) = \sum_{i=1 }^{n}\frac{\ol{u_i}}{i} x^i.
\]
The entries of $A_n$ in $A_n U_n = F_n$ are then
\[
a_{ij} = a(\varphi_i,\varphi_j) = \frac{1}{i+j-1}.
\]
The matrix $A_n$ is called a Hilbert matrix, and has extremely bad condition. For example, for $n=10$, we have $\textrm{cond}(A)=10^{13}$.

Our demands on the choice of the basis functions are \newline
\begin{itemize}
    \item The computation of the matrix elements $a_{ij}=a(\varphi_i,\varphi_j)$ and the right hand side $f(\varphi_j)$ should be "cheap".
    \item In general, we will have $n \gg 100$ and hence $A$ should be sparse.
    \item The matrix $A$ should not be bad conditioned, i.e. $\textrm{cond}(A)\approx O(n)$ or $O(n^2)$ would be okay, but $\textrm{cond}(A)\approx O(n^4)$ or $O(e^n)$ would not be okay.
\end{itemize}

We take $V_n = \spann \l\{ \varphi_i \r\}_{i=1}^n $ as the space of the linear functions
\[
\varphi_j(x) = \begin{cases} \frac{x-x_{j-1}}{h} & \textrm{ if } x\in (x_{j-1},x_j]  \\ \frac{x_{j+1}-x}{h} & \textrm{ if } x \in (x_j,x_{j+1}) \\ 0 & \textrm{ otherwise }\end{cases}.
\]
Here $\l\{ x_j \r\}_{j=0}^n$ is an equidistant mesh of $(0,1)$ with meshsize $h$, i.e. $x_j = j \cdot h$ for $j=0,\dots, N$ with $h=\tfrac{1} N$.\footnote{Notice, that we have $\varphi_i(x_j) = \delta_{ij}$. This property usually is called the \emph{Kronecker-$\delta$-property}.} It then follows, that
\[
a(\varphi_i,\varphi_j) = \begin{cases} \frac{2}{h} & \textrm{ if } i=j \\ -\frac{1}{h} & \textrm{ if } |i-j| = 1 \\ 0 & \textrm{ otherwise } \end{cases},
\]
so $A_n$ is a tridiagonal matrix.

The basis functions hence should have a small support, since this has an impact on the bandwidth of the matrix $A_n$.

\subsubsection*{Design of Finite Element spaces}

\begin{definition}
Let $\Omega \subseteq \R^2$ have a polygonal boundary. A set $\mathcal{T}$ of triangles is called a \textbf{triangulation of $\Omega$} if it holds
\begin{enumerate}[(i)]
    \item
    \[
    \bigcup_{t \in \mathcal{T}} = \ol{\Omega}
    \]
    and
    \item the intersection of two triangles is either a common edge, a point, or empty.
\end{enumerate}
\end{definition}

\begin{example}
{\rot \texttt{The drawings are missing.}}
\end{example}

\begin{definition}
We define
\[
S^{(m)}:= \l\{ v\in C(\ol\Omega) ~|~ v|_t \in \P_m  ~ \forall t \in \mathcal{T} \r\}
\]
as the \textbf{space of the (polynomial) finite elements of order $m$}.
\end{definition}

\begin{remark}\hfill
\begin{itemize}
    \item In general, we have $S^{(m)}\not\subseteq C^1(\Omega)$.
    \item In the case of $m=1,2,3$ we call the finite elements \emph{linear},\emph{quadratic} or \emph{cubic} finite elements.
\end{itemize}
\end{remark}

It is unclear how to construct a basis for $S^{(m)}$. {\rot \texttt{The drawings for $m=1,2,3$ are missing.}} The following lemma tells us how to define basis functions.

\begin{lemma}
Let $m\geq 0$. In a triangle $t$ we have on $m+1$ parallel lines
\[
s= 1+2+\dots+(m+1) = \frac{(m+1)(m+2)}{2}
\]
nodes $\mathcal{N}^{(m)}(t) = \l\{ z_1,\dots,z_s \r\}$. Then, the interpolation problem
\[
\textrm{"Find } p \in \P_m \textrm{ such that } p(z_i)=p_i \textrm{ holds for all }i=1,\dots,s \textrm{."}
\]
has a unique solution for all node values $p_i$.
\end{lemma}
\begin{proof}
The proof will be given later in the lecture.
\end{proof}



Putting the nodes of all triangles together gives us the following corollary.
\begin{corollary}
Let
\[
\mathcal{N}^{(m)} := \bigcup_{t\in \mathcal{T}} \mathcal{N}^{(m)}(t)
\]
be  the collection of all nodes. Then, the interpolation problem
\[
\textrm{"Find } v \in S^{(m)} \textrm{ such that } v(p) = v_p \textrm{ holds for all } p \in \mathcal{N}^{(m)} \textrm{."}
\]
has a unique solution for all node values $v_p$.
\end{corollary}

\begin{definition}
Solving
\[
\textrm{"Find }\varphi_p \in S^{(m)} \textrm{ such that } \varphi_p(q) = \delta_{pq} \textrm{ holds for all } q \in \mathcal{N}^{(m)} \textrm{."}
\]
for each $p\in \mathcal{N}^{(m)}$ gives us the \textbf{nodal basis}
\[
\Lambda^{(m)} := \l\{ \varphi_p ~|~ p \in \mathcal{N}^{(m)} \r\}
\]
of $S^{(m)}$. If $p$ is a vertex of $t$, then $\varphi_p|_t \in \P_m$ is called a \textbf{form function}. We also define
\[
S_h:= \l\{ v \in S^{(m)}~|~ v|_{\partial \Omega} = 0 \r\},
\]
where
\[
h:= \max_{t \in \mathcal{T}} \textrm{diam}(t).
\]
\end{definition}

The next theorem gives a characterization of the space $S_h$.
\begin{theorem}
The space $S_h$ is a closed subspace of $H_0^1(\Omega)$.
\end{theorem}
\begin{proof}
No proof.
\end{proof}

This theorem is good news, since it tells us that $S_h$ is the finite dimensional subspace we were looking for earlier, i.e. we can choose $V_n:= S_h$.

{\rot \texttt{Slides are missing.}}

Notice, that $h$ refers to the size of the triangles, the number of basis functions and maybe also the approximation quality. To precisely measure the approximation quality, we need error estimates. \hfill
\begin{itemize}
    \item Discretization error: $\norm{u-u_h}_{H^1(\Omega)}$, where $u$ is the true solution.
\end{itemize}
If $H_0^1(\Omega)$ is replaced by $S_h$, then we know by Cea's lemma that we only have to consider the \emph{approximation error}:
\begin{itemize}
    \item Approximation Error: $\inf_{v\in S_h} \norm{u-v}_{H^1(\Omega)}$.
\end{itemize}

We now consider the easiest case, i.e. $m=1$ (linear polynomials), $\Omega=(a,b)\subseteq \R$ and
\[
\mathcal{T}: ~ a=x_0 < x_1 < \dots < x_N = b, ~ t_i:= [x_{i-1},x_i], ~ h_i:= x_i - x_{i-1}, ~ h:= \max_{1\leq i \leq N} h_i
\]
where $S_h$ is the space of the continuous and piecewise linear functions with zero bounding conditions.

\begin{theorem}
For each $u\in H_0^1(\Omega)$ it holds
\[
\inf_{v\in S_h} \norm{u-v}_{H^1(\Omega)} \overset{h \to 0}{\longrightarrow} 0
\]
\end{theorem}
\begin{proof}
No proof.
\end{proof}

But now, we are interested in how fast the approximations approach $u$.

We now assume $u\in H^2(a,b)$. By Sobolev's Embedding theory it holds
\[
H^2(a,b) \subseteq C[a,b].
\]
This also holds if $\Omega\subseteq \R^2$ or $\Omega \subseteq \R^3$, but not for $\Omega \in \R^d$ for $d\geq 4$.

Now denote $\norm{\cdot}_{1,\Omega}:=\norm{\cdot}_{H^1(\Omega)}$ and consider
\[
u''(x) = f(x),~ x \in (a,b)=:\Omega, ~ u(a)=u(b)=0.
\]
We now want to compute the discretization error, which we reduce to the estimatia of the interpolation error.
\begin{enumerate}[\textrm{Step }1.]
    \item It holds
    \[
    \norm{u-I_h u}_{1,\Omega}^2 = \sum_{i=1}^{N} \norm{u-I_h u}_{1,t_i}^2.
    \]
    \item We transform onto the unit interval. For each fixed $i$ we transform the unit interval $T=[0,1]$ affine onto $t_i$, i.e.
    \begin{align*}
    t_i & = F_i(T), \\
    x & = F_i(\xi) = x_{i-1} + h_i \xi \\
    \xi_i & := F^{-1}(x) = h^{-1}(x-x_i)
    \end{align*}
    and we set
    \[
    \hat v(\xi) := v(F_i(\xi)) = v(x)
    \]
    for all $v\in H^2(t_i)$.
    {\rot \ttfamily The drawing illustrating the transformation is missing.} By the transformation rule for integrals, it holds\footnote{In the literature, it is common to denote \[\norm{\cdot}_{L^2} = \norm{\cdot}_{H^0} = \norm{\cdot}_{0}.\]}
    \[
    \norm{v}_{0,t_i}^2 = \int_{x_{i-1}}^{x_i}v(x)^2 ~ dx = h_i \cdot \int_{0}^{1}\hat v (\xi)^2 ~ d\xi = h_i \norm{\hat v}_{0,T}^2.
    \]
    The chain rule gives\footnote{Note that we use the chain rule for the strong derivative and transfer it back to the weak derivative by a density argument.}
    \[
    \hat{v}'(\xi) = \frac{d \hat v}{d\xi} (\xi) = \frac{dv}{dx}(x) \cdot \frac{dx}{d\xi} = h_i v'(x)
    \]
    and thus
    \[
    \norm{v'}_{0,t_i}^2 = h_i^{-1}\norm{\hat v'}^2_{0,T}.
    \]
    For $h_i\leq 1$ it holds
    \[
    \norm{v}^2_{1,t_i} \leq h_i^{-1} \norm{\hat v}^2_{1,T} \labtag{III.8}
    \]
    for all $v\in H^2(t_i)$.
    \item We compute the local interpolation error. Applying \eqref{III.8} onto $u- I_h u$ gives
    \begin{align*}
    \norm{u-I_h}^2_{1,t_i} \leq h^{-1}\norm{\hat u - \hat I \hat u_i}^2_{1,T},
    \end{align*}
    where
    \[
    \hat I v(x) = v(0) + (v(1)-v(0))\cdot x
    \]
    is the interpolation operator on $T$. It can then be shown, that
    \[
    \norm{v-\hat I v}^2_{1,T} \leq c \norm{v''}^2_{0,T} \labtag{III.9}
    \]
    with $c= \tfrac 1 3 (8+2 \sqrt{3})$ holds for all $v\in H^2(T)$.
    \item We do the backtransformation. Analogue to step 2, it can be shown by the chain rule, that
    \[
    \norm{\hat u''}_{0,T}^2 = h_i^3 \norm{u''}_{0,t_i}^2.
    \]
    Summing up, we have
    \begin{align*}
    \norm{u - I_h u}^2_{1,\Omega}
    & = \sum_{i=1}^{N}\norm{u- I_h u}^2_{1,t_i} \\
    & \leq \sum_{i=1}^{N}h_i^{-1}\norm{\hat u_i - \hat I \hat u_i}^2_{1,t_i} \\
    & \leq c \cdot \sum_{i=1}^{N} h_i^{-1}\norm{u_i''}^2_{0,T} \\
    & = c \cdot \sum_{i=1}^{N} h_i^{-1} h_i^3 \norm{u''}_{0,t_i}^2 \\
    & \leq c h^2 \norm{u''}_{0,\Omega}^2.
    \end{align*}
    We have thus shown, that for $u\in H^2(a,b)$ the a priori error estimate is given by
    \[
    \norm{u- u_h}_1 \leq \tilde c h \norm{u''}_0 \leq \tilde c h \norm{u}_2
    \]
    with $\tilde c = \sqrt{c} \cdot \tfrac{\Gamma}{\gamma}$.
\end{enumerate}

\begin{remark}[FE in higher dimensions, i.e. $d=2,3$]
{\rot This remark is missing, since it heavily relies on illustrations. Just note, that it can happen that we choose bad triangles (the angles are important, "angle condition") for which the transformation gets bad. The importance of the angles marks a meeting point between geometry and numerical analysis.}
Summing up, the quality of the grid influences the approximation quality.
\end{remark}

\begin{remark}
We are using triangles, since the transformation $F$ gets harder when using different shapes (e.g. squares) since one has to be careful not to flip points.
\end{remark}

\setcounter{subsection}{2}
\subsection{Finite Differences}
Our idea is to replace the derivative of the differential equation by suitable differences.

Consider the interval $[0,1]$ and $x_i=i\cdot h,~ h:= \tfrac 1 N,~ i=0,\dots, N$ and the grid $w_h:= \l\{ x_i ~|~ i=0,\dots, N \r\}$.

\begin{definition}[Gridfunction]
A vector $u_h = \begin{bmatrix} u_0 & \dots & u_N \end{bmatrix}^T \in \R^{N+1}$ assigning each grid point a function value is called \textbf{grid function}.

The restriction $R_h u$ of a function $u\in C^1[0,1]$ onto a grid function is given as
\[
R_h u = \begin{bmatrix} u(x_0) & \dots & u(x_N) \end{bmatrix}^T.
\]
\end{definition}

\begin{example}
Consider the grid $\l\{0.00, 0.25, 0.50, 0.75, 1.00 \r\}$, $u(x) = x^2$ and $R_h u = \begin{bmatrix} 0 & \tfrac 1 {16} & \tfrac 1 4 & \tfrac 9 {16} & 1 \end{bmatrix}^T$.

Caution: Consider $u(x) = \sin(4\pi x)$ and $v(x) = 0$. Then
\[
R_h u = R_h v = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \end{bmatrix}^T.
\]
This means that different functions can have the same grid function.  In this example, we would need a finer mesh for $u$.
\end{example}

\begin{definition}
Let $v(x)$ be a sufficiently smooth function. Denote $v_i := v(x_i)$, where $x_i$ is a node of the grid. Then, we call
\begin{enumerate}[(i)]
    \item
    \[
    D^+ v(x_i) = \frac{v_{i+1}-v_i}{h}
    \]
    the \textbf{forward difference},
    \item
    \[
    D^- v(x_i) = \frac{v_i - v_{i-1}}{h}
    \]
    the \textbf{backward difference},
    \item
    \[
    D^0 v(x_i) = \frac{v_{i+1}-v_{i-1}}{2h}
    \]
    the \textbf{central difference} and
    \item
    \[
    D^+ D^- v(x_i) = \frac{v_{i+1}-2v_i + v_{i-1}}{h^2}
    \]
    the \textbf{second difference}.
\end{enumerate}
\end{definition}

\begin{definition}
Let $L$ be a differntial operator. Then the difference operator $L_h:\R^{N+1}\to \R^N$ is \textbf{consistent with $L$ of order $k$} if
\[
\max_{0\leq i \leq N} \norm{(Lu)(x_i) - (L_h u_h)_i} =: \norm{(Lu)(x_i) - (L_h u_h)_i}_{\infty,d} = O(h^k)
\]
holds.
We call $\norm{\cdot}_{\infty,d}$ the \textbf{discrete maximum norm in the space of the grid functions}.
\end{definition}

\begin{remark}
Consistency is a measure for the approximation quality.
\end{remark}

\begin{example}
The operators $D^+, D^-,D^0$ are consistent of order $1$ to $L=\tfrac{d}{dx}$ and $D^+D^-$ is consistent of order $2$ to $L=\tfrac{d^2}{dx^2}$.
\end{example}

\begin{example}
Consider
\[
Lu = \frac{d}{dx} \l( k(x)\cdot \frac{du}{dx} \r),
\]
where $k(x)$ is continuously differentiable. We define $L_h$ by
\begin{align*}
    (L_h u_h)_i
    & = D^+ (a D^- u(x_i)) \\
    & = \frac{1}{h}\l( a_{i+1} D^- u (x_{i+1}) - a_i D^- u(x_i)\r) \\
    & = \frac{1}{h} \l( a_{i+1} \frac{u_{i+1}-u_i}{h}- a_i \frac{u_i - u_{i-1}}{h}\r),
\end{align*}
where $a$ is a grid function, that has to be chosen suitable (cp $k(x)$). Using the chain rule, we obtain
\[
(Lu)_i = k'(x_i) (u')_i + k(x_i) + (u'')_i.
\]
A Taylor expansion gives
\[
(L_h u_h)_i = \frac{a_{i+1}-a_i }{h}(u')_i + \frac{a_{i+1} - a_i}{2}(u'')_i + \frac{h (a_{i+1} - a_i)}{6} (u''')_i + O(h^2).
\]
For the difference we obtain
\begin{align*}
&(Lu)_i - (L_h u_h)_i \\ = &\l(k(x_i) - \frac{a_{i+1}-a_i }{h}\r) (u')_i + \l( k(x_i) - \frac{a_{i+1} - a_i}{2} \r) (u'')_i  \\  & - \frac{h (a_{i+1}-a_i)}{6} (u''')_i + O(h^2).
\end{align*}
$L_h$ is consistent with $L$ of order $2$, if
\[
\frac{a_{i+1}-a_i}{h} = k'(x_i) + O(h^2) ~\und~ \frac{a_{i+1}-a_i} 2 = k(x_i)+ O(h^2).
\]
Possible choises are
\[
a_i = \frac{k(x_i)-k(x_{i+1})}{2}, ~ a_i = k \l( x_i - \tfrac{h}2 \r) ~\und~ a_i = \sqrt{k(x_i)\cdot k(x_{i-1})}.
\]
Note, that the choice $a_i=k(x_i)$ gives only order of consistency $1$.
\end{example}

Now consider
\[
Lu = u'' + b(x) u' + c(x) u = f(x), ~ x \in (0,1) \labtag{III.10}, ~ u(0)=u(1)=0
\]
where we write $u=u(x)$ for short. We assume, that $b$ and $c$ are sufficiently smooth and that $c(x)\geq 0$ for all $x\in[0,1]$.

\begin{definition}
The central difference scheme for \eqref{III.10} is given as
\[
- D^+ D^- u_i + b_i D^0 u_i + c_i u_i = f_i, ~ i=1,\dots, N, ~ u_0 = u_N = 0.
\]
\end{definition}

\begin{remark}
This leads to a tridiagonal matrix system of linear equations
\[
r_i u_{i-1} + s_i u_i + t_i u_{i+1} = f_i, ~ i=1\dots,N, ~ u_0 = u_N = 0.
\]
with
\[
r_i = - \frac{1}{h^2} - \frac{1}{2h}b_i, ~ s_i = c_i + \frac{2}{h^2} \und t_i = -\frac{1}{h^2} + \frac{1}{2h}b_i.
\]
\end{remark}

Transferring the consistency of Definition III.21 to \eqref{III.10}, we obtain
\[
L_h u_h := R_h (Lu).
\]
We assume, that the boundary conditions are imposed such that the first row and the last row of $L_h$ are $\begin{bmatrix} 1 & 0 & \dots & 0 \end{bmatrix}^T$ and $\begin{bmatrix} 0 & \dots & 0 & 1 \end{bmatrix}^T$, respectively. It holds
\[
\l( R_h(Lu) \r)_0 = u_0, ~ \l( R_h(Lu) \r)_N = u_N.
\]
We thus call \eqref{III.10} consistent of order $k$, if
\[
\norm{L_h R_h u - R_h (Lu)}_{\infty,d} \leq c\cdot h^k,
\]
where $c$ and $k$ are independent of $h$.\footnote{For $u \in C^4[0,1]$ it can be shown, that the difference scheme for \eqref{III.10} has order $k=2$.}

\begin{definition}
A difference scheme $L_h u_h = f_h$ is called \textbf{stable}, if there exists a stability constant $c_s$, which is not depending on the meshsize $h$, such that
\[
\norm{u_h}_{\infty,d} \leq c_s \norm{L_h u_h}_{\infty,d}
\]
holds for all grid functions $u_h$.
\end{definition}

\begin{definition}
A difference scheme for \eqref{III.10} is \textbf{convergent of order $k$} if there exist constants $c,k>0$ (being independent of $h$) such that
\[
\norm{u_h - R_h u}_{\infty,d} \leq c\cdot h^k.
\]
\end{definition}

For a consistent and stable scheme we have
\begin{align*}
\norm{u_h - R_h u}_{\infty,d}
& \leq c_s \norm{L_h(u_h - R_h u)}_{\infty,d} \\
& = c_s \norm{L_h u_h - L_h R_h u}_{\infty,d} \\
& = c_s \norm{f_h - L_h R_h u}_{\infty,d} \\
& = c_s \norm{R_h f - L_h R_h u}_{\infty,d} \\
& = c_s \norm{R_h L u - L_h R_h u}_{\infty,d} \\
& \leq K \cdot h^k.
\end{align*}

We have thus shown the following theorem.

\begin{theorem}
A consistent and stable difference scheme is convergent and the order of consistency and convergence are the same.
\end{theorem}

\begin{remark}
The consistency is usually shown by Taylor expansion. The stability is in general harder to show, dealing with functions and matrices.
\end{remark}

\begin{definition}
Let $x,y \in \R^n$ and $A=[a_{ij}]\in\R^{n,n}$. Then we define
\begin{align*}
x \leq y &~~:\Leftrightarrow ~~x_i \leq y_i ~ \forall i=1,\dots,n \\
x \geq 1 &~~:\Leftrightarrow ~~x_i \geq 1 ~ \forall i=1,\dots,n \\
A \geq 0 &~~:\Leftrightarrow ~~a_{ij} \geq 0 ~ \forall i,j = 1,\dots,n.
\end{align*}
We call $A$ \textbf{inverse monotone}, if $A^{-1}$ exists and $A^{-1}\geq 0$.
\end{definition}

Let us assume, that $A$ is inverse monotone and $Av \leq Aw$ for $v,w \in \R^n$. Then
\[
A(v-w) = b \leq 0
\]
and thus
\[
v-w = A^{-1} b \leq 0
\]
so $v\leq w$. \newline \hfill

Now, consider the BVP
\[
\begin{aligned}
(Lu)(x) = & - (pu')'(x) + q(x) u'(x) + r(x)u(x) = f(x), ~~ x \in I := [a,b] \\
& u(a) = \alpha, ~ u(b) = \beta,
\end{aligned} \labtag{III.11}
\]
the point grid
\[
a = x_0 < x_1 < \dots < x_N = b, ~ I_i=[t_{i-1},t_i)
\]
and the difference scheme
\begin{align*}
L_h u_i := & - D_{h/2}^0 (p_i D^0_{h/2}u_i) + q_i \cdot D_h^0 u_i + r_i u_i = f_i, ~ i=1,\dots,N-1 \\
& u_0 = \alpha, ~ u_N = \beta,
\end{align*}
where
\[
D^0_h u(x) = \frac{u(x+h)-u(x-h)}{2h} \und g_i = g(x_i).
\]
This is equivalent to
\[
\tilde A_h \tilde u_h = \tilde f_h \labtag{III.12}
\]
where $\tilde A_h$ and $\tilde f_h$ are the result of $u_0 = \alpha, u_N = \beta$ and
\[
- \l( p_{i-1/2} u_i - (p_{i-1/2}+p_{i+1/2})y_i + p_{i+1/2}u_{i+1} + \tfrac h 2 q_i (u_{i+1}-u_{u-1}) + h^2 r_i y_i \r) = h^2 f_i.
\]
{\rot $A_h u_h = f_h$}, where
\[
A_h = \einsdurch{h^2} \begin{bmatrix}
p_{1/2} + p_{3/2} + h^2 r_1 & - p_{3/2} + \tfrac{1}2 h q_1 & &  \\
- p_{3/2} + \tfrac 1 2 h q_N & \ddots & \ddots & \\
& \ddots & \ddots & \ddots \\
\end{bmatrix}
\]
and
\[
f_h = \l( f_1 + h^{-2} p_{1/2}\alpha + \tfrac 1 2 h^{-1} q \alpha, ~ f_2, \dots, f_{N-1} + h^{-1} p_{N-1/2} \beta - \tfrac 1 2 h^{-1} q_N \beta \r),
\]
where $A_h$ is triadiagonal.

\begin{remark}\hfill
\begin{itemize}
    \item It can be shown that this scheme for \eqref{III.11} has order of consistency $2$.
    \item Stability in difference schemes is in general harder to show. {\rot The similarity of $A_h$ to $L_h$ is missing.}
\end{itemize}
\end{remark}

\begin{definition}
A matrix $A=[a_{ij}]\in\R^{N\times N}$ is called a \textbf{M-matrix} if\hfill
\begin{enumerate}[(a)]
    \item $a_{ij}\leq 0$ for $i\neq j$
    \item $A^{-1}$ exists with $A^{-1}\geq 0$.
\end{enumerate}
\end{definition}

Let us assume, that in \eqref{III.11} we have $p>0, q = 0$ and $r\geq 0$. Then the matrix $A_h$ is given by
\[
A_h = \einsdurch{h^2} \begin{bmatrix}
p_{1/2} + p_{3/2} + h^2 r_1 & - p_{3/2} & & \\
-p_{i-1/2} & \ddots & \ddots & \\
& \ddots & \ddots & \ddots
\end{bmatrix}.
\]
This matrix is then \emph{strictly diagonal dominant}, i.e. it holds for at least one $s\in \l\{ 1, \dots, N \r\}$ that
\[
\sum_{j\neq i} |a_{ij}| \leq |a_{ii}|, ~ 1 \leq i \leq N, ~\und~ \sum_{j\neq s} |a_{sj}| < |a_{ss}|.
\]
Further, it is \emph{irreducable}, i.e. for each two indices $i,j \in \l\{ 1,\dots, N \r\}$ there exists a sequence of indices $j_1,\dots, j_n$ such that
\[
a_{j_1,i} \neq 0, ~ a_{j_2,j_1} \neq 0, \dots, ~ a_{j,j_m} \neq 0.
\]

\begin{theorem}
It can be shown, that for \eqref{III.11} with $q=0,r \geq 0$ the difference scheme \eqref{III.12} is stable.
\end{theorem}

\footnote{If interested, the books \emph{R. Plato, "Numerische Mathematik - kompakt"} and \emph{St\"{o}r/Bullirsch, "Numerical Mathematics II"} can give further information and proofs regarding difference methods.}

In case of $q\neq 0$, we have
\[
A_h = \einsdurch{h^2} \begin{bmatrix}
p_{1/2} + p_{3/2} + h^2 r_1 & -p_{3/2} + 1/2 h q_1 & & \\
- p_{i-1/2} - 1/2 h q_i & \ddots & \ddots & \\
& \ddots & \ddots & \ddots
\end{bmatrix}.
\]
This matrix is unsymmetric and only diagonal dominant, if
\[
h \leq 2 \min_{1\leq i \leq N-1} \l\{ \frac{\min \l\{ p_{i-1/2},p_{i+1/2} \r\}}{|q_i|} \r\}.
\]

We now consider the problematic case $|q_i| \gg |p_i|$.

\begin{example}
Let $I=[0,1], ~ f = 0, ~ q = 1, ~ r = 0, ~ 0 < p = \varepsilon \ll 1$ and consider the problem
\begin{align*}
L^{\varepsilon}u(x) = &  - \varepsilon u''(x) + u'(x) = 0, ~ x \in I \\
 & u(0) = 1, u(1) = 0.
\end{align*}
This problem has the analytical unique solution
\[
u^{\varepsilon}(x) = \frac{e^{1/\varepsilon} - e^{x/\varepsilon}}{e^{1/\varepsilon}-1}.
\]
For $x = 1 - \delta, ~ \delta > 0$, we have
\[
u^\varepsilon(1-\delta) = \frac{e^{1/\varepsilon}}{e^{1/\varepsilon}-1} \l( 1 - e^{-\delta/\varepsilon} \r) \approx 1.
\]
Further, for $\varepsilon=0$ we have $u_0 = 1$, which does \textbf{not} meet the boundary condition for $x=1$.

{\rot The drawing illustrating the problem with the solution is missing.}

Approximating this problem by the above difference scheme with equidistant step size $h=\tfrac 1 N$ we have
\[
\begin{aligned}
    - (\varepsilon + \tfrac 1 2 h) u_{i-1} + 2 \varepsilon u_i - (\varepsilon - \tfrac 1 2 h) u_{i+1}  & = 0, ~ 1 \leq i \leq N-1 \\
    u_0 = 1, & u_N = 0
\end{aligned} \labtag{III.13}
\]
The corresponding system matrix is only non-negative and diagonal dominant, if $h\leq 2 \varepsilon$.

For the solution of the difference equation \eqref{III.13} we now take the approach $u_i = \lambda^i$ (cp. with characteristic polynomials of chapter II). We are seeking for the roots of the system. It holds
\[
\lambda^i + \frac{2 \varepsilon}{\tfrac h 2 \varepsilon} \lambda - \frac{\tfrac h 2 + \varepsilon}{\tfrac h 2 - \varepsilon} = 0, ~ u_0 = 1, u_N = 0.
\]
We take the ansatz
\[
u_i = c_1 \lambda_1^i + c_2 \lambda_2^i,
\]
where $\lambda_1,\lambda_2$ are the roots of the polynomial.
It can then be shown, that the solution is given as
\[
u_i = \frac{\lambda_1^N \cdot \lambda_2^i - \lambda_2^N \cdot \lambda_1^i}{\lambda_1^N - \lambda_2^N}
\]
and the roots are
\[
\lambda_1 = 1 \und \lambda_2 = \frac{\varepsilon + \tfrac h 2}{\varepsilon - \tfrac h 2},
\]
so for $\varepsilon \ll \tfrac h 2$ we have $\lambda_2 \approx -1$. In this case, we would have an oscillating solution
\[
u_i = \frac{\lambda_2^i - \lambda_2^N}{1-\lambda_2^N}
\]

We can get rid of this, by \hfill
\begin{enumerate}[(1).]
    \item Compute $u'(t)$ by $D_h^+$ or $D_h^-$. Because for
    \[
    D_h^- u(t) = \frac{u(t+h) - u(t)}{h}
    \]
    we arrive at
    \[
    -(\varepsilon + h) u_{u-1} + (2 \varepsilon + h) u_i - \varepsilon u_{i+1} = 0.
    \]
    The corresponding system matrix is diagonal dominant for all $h>0$ and is a M-matrix. With $u_i = \lambda^i$ we arrive at the polynomial
    \[
    \lambda^2 - \frac{2 \varepsilon + h}{ \varepsilon} \lambda + \frac{\varepsilon+h}{\varepsilon} = 0,
    \]
    which has the roots
    \[
    \lambda_1 = \frac{\varepsilon + h}{\varepsilon} ~\und~ \lambda_2 = 1,
    \]
    so the solution is
    \[
    u_i = \frac{\lambda_1^N - \lambda_1^i}{\lambda_1^N - 1}.
    \]
    This is called the $D_h^-$ upwind discretization (\emph{R\"{u}ckw\"{a}rtsdiskretisierung}).
    \item Artificial diffusion. We set $\hat \varepsilon := \varepsilon + \delta h$ and obtain
    \[
    -(\hat \varepsilon + \tfrac h 2) u_{i-1} + 2 \hat \varepsilon u_i - (\hat \varepsilon - \tfrac h 2) u_{i+1} = 0,
    \]
    which works for $\delta \geq 1/2$, $A_h$ is a M-matrix.
\end{enumerate}
\end{example}

Reconsider the problem \eqref{III.11}. If we have $p=1$ and $q=r=0$, it follows that
\[
A_h = \einsdurch{h^2} \begin{bmatrix}
2 & -1 & & \\
- 1 & \ddots & \ddots & \\
& \ddots & \ddots & \ddots
\end{bmatrix}.
\]
Since we have
\[
\textrm{cond}(A_h) = \frac{\lambda_{\textrm{max}}(A_h)}{\lambda_{\textrm{min}}(A_h)}
\]
and the eigenvalues and eigenvectors are
\[
\lambda_k = \einsdurch{h^2} \l\{ 2 - 2 \cos(kh\pi) \r\}, ~ k=1,\dots,N
\]
and
\[
w^k = \l\{ \sin{ikh\pi} \r\}_{i=1}^N,
\]
respectively, it holds
\[
\lambda_{\textrm{max}} = \frac{1}{h^2} (2 - 2 \cos((1-h)\pi)) = \frac{4}{k^2} + O(1)
\]
and
\[
\lambda_{\textrm{min}} = \frac{1}{h^2} (2 - 2 \cos(h\pi)) = \frac{1}{h^2} (2-2(1- \tfrac 1 2 h^2 \pi^2)) + O(h^4).
\]
Thus, we have
\[
\textrm{cond}(A_h) = \frac{\lambda_{\textrm{max}}}{\lambda_{\textrm{min}}} = \frac{4}{h^2 \pi^2}+ O(1).
\]
Hence, the condition is $O(h^{-2})$ and thus the finer the discretization, the worse is the condition of $A_h$.

\subsubsection*{Box-Method}
There exist also so called Box methods for finite difference methods. Here, the node values $v_i$ of $u$ are representable for the value of $u$ in an interval, i.e. $[x_{i-1/2},x_{i-1/2}] = [x_i - \tfrac{h}{2}, x_i + \tfrac{h}2]$ and we assume that $u_i$ is piecewise constant on this interval.

\begin{example}
Consider
\[
- (k(x)u')' + c(x) u(x) = f(x), ~ 0 < x < 1, ~ u(0) = u(1) = 0, ~ k(x),c(x) \geq 0.
\]
Then we want that
\begin{align*}
\int_{x_{i-1/2}}^{x_{i+1/2}}-(k(x) u')' + c(x) u(x) ~ dx
& = - (k(x) u'(x))\bigg|^{x_{i+1/2}}_{x_{i-1/2}} + \int_{x_{i-1/2}}^{x_{i+1/2}} c(x) u(x)~ dx \\
& \overset{\textrm{!}}{=} \int_{x_{i-1/2}}^{x_{i+1/2}}f(x)~ dx
\end{align*}
holds. Hence, we can derive
\[
\einsdurch{h} \l(k_{i+1/2} \frac{u_{i+1} u_i}{h} - k_{i-1/2} \frac{u_i - u_{i-1}}{h}\r) + c_i u_i = f_i.
\]
\end{example}

\subsubsection*{Shooting Methods}
Given some BVP
\begin{align*}
y' & = f(x,y) \in \R^n, \\
a & \leq x \leq b, \\
R(y(a),y(b)) & = 0 \in \R^n,
\end{align*}
we take some guess $s:=y(a)$ and solve the IVP
\[
y' = f(x,y), ~ y(a) = s.
\]
The solution $y(x;a,s)$ of this IVP is then evaluated at the right hand side of the boundary
\[
F(s):= R \Big( y(a;a,s),~y(b;a,s) \Big) = R \Big( s, y(b;a,s) \Big).
\]
With regard to the BVP, we seek for a root of $F$.

{\rot The drawing explaining shooting methods is missing.}

We do single shooting by using an inexact Newton method. We are seeking for the roots of
\begin{align*}
    DF(s^{(i)}) \Delta s^{(i)} & = - F(s^{(i)}), \\
    s^{(i+1)} & = s^{(i)} + \lambda^{(i)} \Delta s^{(i)}, \\
    \Delta s^{(i)} & = s^{(i)} - s^{(i-1)},
\end{align*}
where $DF$ is a Jacobian, which is not computed exactly, but numerical.

We are using following steps:

\begin{enumerate}[(1).]
    \item[(0).] Guess $s^{(i)},~i=0,~\lambda^{(-1)}:= 1$ and compute $F(s^{(i)})$.
    \item Compute $DF s^{(i)}$.
    \item Compute $\Delta s^{(i)},~ \lambda^{(i)} = \min \l\{ 1,2 \lambda^{(i-1)} \r\}$.
    \item Set $s^{(i+1)} := s^{(i)} + \lambda^{(i)} \Delta s^{(i)}$.
    \item Compute $F(s^{(i+1)})$.
    \item If $\norm{F(s^{(i+1)})} < \textrm{TOL}$, then end. Otherwise, go to (6).
    \item If $\norm{F(s^{(i+1)})} > \norm{F(s^{(i)})}$, set $\lambda^{(i)} = \tfrac{\lambda^{(i)}} 2$ and go to (3). Otherwise, set $i=i+1$.
\end{enumerate}

\begin{remark}\hfill
\begin{itemize}
    \item Evaluating $F$ is equivalent to solving an IVP.
    \item We have
    \[
    DF = \partiell{R\Big( y(a;a,s),~y(b;a,s) \Big)}{y_a} + \partiell{R \Big( y(a;a,s), ~ y(b;a,s) \Big)}{y_b} \cdot \underbrace{\partiell{y(b;a,s)}{s}}_{=: G}.
    \]
    The matrix $G$ is called \emph{sensitivity matrix}.
\end{itemize}
\end{remark}

The sensitivity matrix $G$ can be approximated via a difference approximation. More precisely, the $i$-th column of $G$ is given by
\[
\partiell{y(x;a,s)}{s_i} \approx \frac{y(x;a,s+\delta e_i)-y(x;a,s)}{\delta}, \tag{$\ast$}
\]
where $\delta$ is a discretization parameter.

It is crucial to choose $\delta$ properly. However, it is unclear how to do that.

We have two kinds of errors: \hfill
\begin{itemize}
    \item the approximation error, which is $O(\delta)$ and
    \item the rounding error, since the smaller $\delta$ is, the closer are $y(x;a,s+\delta e_i)$ and $y(x;a,s)$.
\end{itemize}
In case of a linear BVP, the Newton method converges after $1$ step. However, in the case of a nonlinear BVP, it can happen that for a bad guess $s$ the Newton scheme does not converge, even if the BVP has locally a unique solution.

\begin{example}
Consider
\[
y'=y^2, ~ y=\frac 1 {c-t}, ~ c\in\R
\]
with the boundary condition $y(10) -y(0) = \tfrac{1}{20}$. Then $c=20,~y(t) = \tfrac{1}{20-t}$. It follows that $s = y(0) = \tfrac{1}{20}=0.05$.  Taking a "bad" initial guess, e.g. $s^{(0)}=y^{(0)}(0)=0.2$, we get $y=\tfrac{1}{5-t}$ and thus the corresponding IVP has a singularity in $[0,10]$. Hence, we can not evaluate $F$.

Thus, a good starting guess is important and a single shooting method is not suitable for highly nonlinear problems.
\end{example}



























































































\end{document}