% vim: set sts=4 et :
\documentclass[11 pt,xcolor={dvipsnames,svgnames,x11names,table}]{beamer}
\usetheme[
bullet=circle, % Other option: square
bigpagenumber, % circled page number on lower right
topline=true, % colored bar at the top of the frame
shadow=false, % Shading for beamer blocks
watermark=BG_lower, % png file for the watermark
]{Flip}
\usepackage{bm}
\usepackage{tikz,pgfplots,pgfplotstable}
\usepgfplotslibrary{fillbetween}
\usepackage{graphicx}
% SOME COMMANDS THAT I FIND HANDY
% \renewcommand{\tilde}{\widetilde} % dinky tildes look silly, dosn't work with fontspec
\newcommand{\comment}[1]{\textcolor{comment}{\footnotesize{#1}\normalsize}} % comment mild
\newcommand{\Comment}[1]{\textcolor{Comment}{\footnotesize{#1}\normalsize}} % comment bold
\newcommand{\COMMENT}[1]{\textcolor{COMMENT}{\footnotesize{#1}\normalsize}} % comment crazy bold
\newcommand{\Alert}[1]{\textcolor{Alert}{#1}} % louder alert
\newcommand{\ALERT}[1]{\textcolor{ALERT}{#1}} % loudest alert
%% "\alert" is already a beamer pre-defined
\newcommand*{\Scale}[2][4]{\scalebox{#1}{$#2$}}%
\graphicspath{{images/}} % Put all images in this directory. Avoids clutter.
% suppress frame numbering for backup slides
% you always need the appendix for this!
\newcommand{\backupbegin}{
\newcounter{framenumberappendix}
\setcounter{framenumberappendix}{\value{framenumber}}
}
\newcommand{\backupend}{
\addtocounter{framenumberappendix}{-\value{framenumber}}
\addtocounter{framenumber}{\value{framenumberappendix}}
}
\begin{document}
\tikzstyle{every picture}+=[remember picture]
{
\setbeamertemplate{sidebar right}{\llap{\includegraphics[width=\paperwidth,height=\paperheight]{bubble2}}}
\begin{frame}[c]%{\phantom{title page}}
\begin{center}
\begin{center}
\begin{columns}
\begin{column}{0.9\textwidth}
\flushright%\fontspec{Trebuchet MS}
\bfseries \Huge {Ordinary Differential Equations (I)}
\end{column}
\begin{column}{0.2\textwidth}
%\includegraphics[width=\textwidth]{SHiP-2}
\end{column}
\end{columns}
\end{center}
\quad
\vspace{3em}
\begin{columns}
\begin{column}{0.6\textwidth}
\flushright \vspace{-1.8em} {%\fontspec{Trebuchet MS}
\large Marcin Chrzaszcz, Danny van Dyk\\\vspace{-0.1em}\small \href{mailto:mchrzasz@cern.ch}{mchrzasz@cern.ch}, \href{mailto:danny.van.dyk@gmail.com}{danny.van.dyk@gmail.com}}
\end{column}
\begin{column}{0.4\textwidth}
\includegraphics[height=1.3cm]{uzh-transp}
\end{column}
\end{columns}
\vspace{1em}
\vspace{0.5em}
\textcolor{normal text.fg!50!Comment}{Numerical Methods, \\ 14. November, 2016}
\end{center}
\end{frame}
}
\begin{frame}{Plan for today}
\begin{itemize}
\item \alert{Catching up with last time}\\
\begin{itemize}
\item What does an implicit Runge-Kutta method look like?
\end{itemize}
\vfill
\item \alert{So far we discussed ODEs of first order}\\
\begin{itemize}
\item How to turn an ODE of order $n$ into a system of $n$ ODEs of first order?
\end{itemize}
\vfill
\item \alert{Specific problem of a given ODE in a boundary value problem}\\
\begin{itemize}
\item How can we use our knowledge of initial value problems to solve these?
\item Can we solve these directly, w/o resorting to initial value problems?
\end{itemize}
\end{itemize}
\end{frame}
\begin{frame}{Implicit Runge-Kutta methods}
\begin{itemize}
\item same as before, method of stage $s$ and order $p$
\begin{equation*}
{\color{red}k_i} = f\left[x_n + c_i h, y_n + h \sum_{j} a_{i,j} {\color{red}k_j} \right]
\end{equation*}
\item compare Butcher Tableau, which is now fully populated
\end{itemize}
\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{c|ccccc}
$c_1$ & $a_{1,1}$ & $a_{1,2}$ & $\dots$ & $a_{1,s-1}$ & $a_{1,s}$ \\
$\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ \\
$c_s$ & $a_{s,1}$ & $a_{s,2}$ & $\dots$ & $a_{s,s-1}$ & $a_{s,s}$ \\
\hline
& $b_1$ & $b_2$ & $\dots$ & $b_{s-1}$ & $b_s$ \\
\end{tabular}
\renewcommand{\arraystretch}{1.0}
\end{center}
\begin{itemize}
\item while explicit methods have strictly $s \geq p$, implicit method
can achieve orders $p > s$
\end{itemize}
\end{frame}
\begin{frame}{Gauss-Legendre methods \dots}
\hfill \dots a subclass of implicit Runge-Kutta methods\\
\vfill
Revisit the integral we need:
\begin{equation*}
y_{n + 1} = y_n + {\color{red} \int_{x_n}^{x_n + h} \mathrm{d}x\, f(x, y(x))}
\end{equation*}
\begin{itemize}
\item reuse what we know about Gauss quadratures
\item rescale the problem such that Legendre polynomials are the appropriate
orthogonal basis for $f(x, y(x))$
\item example: Gauss Legendre method of stage $2$:
\end{itemize}
\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{c|cc}
$\frac{1}{2} - \frac{\sqrt{3}}{6}$ & $\frac{1}{4}$ & $\frac{1}{4} - \frac{\sqrt{3}}{6}$ \\
$\frac{1}{2} - \frac{\sqrt{3}}{6}$ & $\frac{1}{4} + \frac{\sqrt{3}}{6}$ & $\frac{1}{4}$ \\
\hline
& $\frac{1}{2}$ & $\frac{1}{2}$ \\
\end{tabular}
\renewcommand{\arraystretch}{1.0}
\end{center}
\begin{itemize}
\item practical aspects of implicit Runge-Kutta methods requires a tangent
\end{itemize}
\end{frame}
\begin{frame}{Tangent: Fixed-point problems}
\alert{What is a fixed-point problem?}\\
\vfill
\begin{itemize}
\item consider an automorphism $f$, e.g.: a function $\mathbb{R}\ni z \mapsto f(z)\in \mathbb{R}$
\item a point $z^*$ is a fixed point under $f$ if it fulfills
\end{itemize}
\begin{equation*}
f(z^*) = z^*
\end{equation*}
\begin{itemize}
\item A \alert{fixed-point problem} is finding the fixed point $z^*$ for a given $f$.
\item Iteratively, we can approach $\lim_{n\to \infty} z_n = z^*$ via:
\end{itemize}
\begin{equation*}
z^{(i + 1)} = f(z^{(i)})
\end{equation*}
\begin{itemize}
\item this works, as long as $z^*$ is an attractive fixed point, i.e.:
\end{itemize}
\begin{equation*}
|| f(z^{(i+1)}) - f(z^{(i)}) || < || z^{(i + 1)} - z^{(i)} ||
\end{equation*}
\begin{itemize}
\item for further reading see the Banach fix-point theorem!
\end{itemize}
\end{frame}
\begin{frame}{Practical Aspects}
We can view implicit single step methods
\begin{equation*}
y' = f(x, y)\qquad y_{n + 1} = y_{n} + \Phi(F; x_n, h; y_n, y_{n+1})
\end{equation*}
as a fixed point problem:
\begin{itemize}
\item start with a guess for $y_{n,i=0}$, e.g. based on an explicit method
\item compute the new approximation of the implict result for $y_{n+1}$ as
\begin{equation*}
y_{n+1, i+1} = y_n + \Phi(F; x_n, h; y_n, y_{n+1, i})
\end{equation*}
\item stop after a fixed number of iterations, or when the difference between
previous and current iteration falls below a pre-determined threshold.
\item check attraction requirement in every step
\end{itemize}
\end{frame}
\begin{frame}{Example}
\begin{itemize}
\item first step: Euler explicit
\item next steps: repeat Euler implicit until $| z^{(i + 1)} - z^{(i)} | < 2\cdot 10^{-3}$
\end{itemize}
\begin{columns}
\begin{column}{.6\textwidth}
\resizebox{\textwidth}{!}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=left,
ymin=-0.40,ymax=+1.30,
xmin=+0.00,xmax=+8.50
]
\addplot[
ultra thick,
gray,
mark=none,
domain=0:8.5
]
{ exp(-0.5) };
\addplot[
thick,
color=ForestGreen,
mark=+,
only marks
]
coordinates {
(1, 0.5)
(2, 0.75)
(3, 0.625)
(4, 0.6875)
(5, 0.65625)
(6, 0.671875)
(7, 0.664063)
(8, 0.667969)
(9, 0.666016)
};
\only<2->{
\addplot[
thick,
color=orange,
mark=+,
only marks
]
coordinates {
(1, 0)
(2, 1)
(3, 0)
(4, 1)
(5, 0)
(6, 1)
(7, 0)
(8, 1)
(9, 0)
};
}
\only<3->{
\addplot[
thick,
color=red,
mark=+,
only marks
]
coordinates {
(1, -0.05)
(2, +1.0525)
(3, -0.105125)
(4, +1.11038)
(5, -0.1659)
(6, +1.1742)
(7, -0.232905)
(8, 1.24455)
(9, -0.306778)
};
}
\end{axis}
\end{tikzpicture}
}
\end{column}
\begin{column}{.4\textwidth}
\begin{itemize}
\item works {\color{ForestGreen}fine} for $f(x, y) = -y$
\item<2-> {\color{orange}fails} for $f(x, y) = -2 y$
\item<3-> {\color{red}fails spectacularly}
for $f(x, y) = -2.1 y$
\end{itemize}
\end{column}
\end{columns}
\only<3->{\alert{Be careful with your starting value!}}
\end{frame}
\begin{frame}{Toward Boundary Problems: ODEs of order $2$}
Let's start with an explicit ODE of order $2$:
\begin{align*}
y'' & = f(x; y, y')
\end{align*}
Let's rename $y(x) \to z_1(x)$ introduce a new function $z_2(x) = y'(x)$, such that
we can write
\begin{align*}
z_2' & = f_2(x; z_1, z_2) & & [\equiv f(x; y = z_1, y' = z_2)]\\
z_1' & = f_1(x; z_1) & & [\equiv z_2]
\end{align*}
We can now solve the (vector-valued) ODE of order $1$.
\end{frame}
\begin{frame}{Generalization to order $N$}
\begin{itemize}
\item we can apply this now to any order $N$ ODE to obtain
a system of $N$ coupled ODEs of order $1$
\item however, a first-order system of $N$ equations
requires $N$ initial conditions to solve!
\item consequently, we already needed $N$ initial conditions
to solve the order $N$ system!
\end{itemize}
\end{frame}
\pgfplotstableread{
x J Jp
1. 0.440051 0.325147
1.05 0.455897 0.308608
1.1 0.470902 0.291529
1.15 0.485041 0.273945
1.2 0.498289 0.255892
1.25 0.510623 0.237407
1.3 0.522023 0.21853
1.35 0.53247 0.199297
1.4 0.541948 0.17975
1.45 0.550441 0.159927
1.5 0.557937 0.13987
1.55 0.564424 0.11962
1.6 0.569896 0.0992172
1.65 0.574344 0.0787044
1.7 0.577765 0.058123
1.75 0.580156 0.0375147
1.8 0.581517 0.0169214
1.85 0.581849 -0.00361515
1.9 0.581157 -0.0240536
1.95 0.579446 -0.0443527
2. 0.576725 -0.0644716
}{\BesselJtable}
\pgfplotstableread{
1. 0.44051 -0.490276
1.05 0.416572 -0.467924
1.1 0.393655 -0.449309
1.15 0.371591 -0.433691
1.2 0.350246 -0.420483
1.25 0.329511 -0.409212
1.3 0.309299 -0.399494
1.35 0.289541 -0.391015
1.4 0.270181 -0.383519
1.45 0.251176 -0.37679
1.5 0.232492 -0.37065
1.55 0.214104 -0.36495
1.6 0.195992 -0.359565
1.65 0.178144 -0.354388
1.7 0.160551 -0.349331
1.75 0.14321 -0.344316
1.8 0.12612 -0.339281
1.85 0.109283 -0.33417
1.9 0.092705 -0.328939
1.95 0.0763921 -0.323548
2. 0.0603534 -0.317967
}{\BesselYtable}
\pgfplotstableread{
1 0.440051 0.325147
1.05 0.455468 0.308347
1.1 0.47002 0.291039
1.15 0.483683 0.273261
1.2 0.496435 0.255049
1.25 0.508257 0.236443
1.3 0.519132 0.21748
1.35 0.529042 0.198202
1.4 0.537974 0.178647
1.45 0.545917 0.158856
1.5 0.55286 0.138869
1.55 0.558797 0.118729
1.6 0.56372 0.0984759
1.65 0.567628 0.078151
1.7 0.570518 0.0577958
1.75 0.57239 0.0374514
1.8 0.573248 0.0171588
1.85 0.573096 -0.00304138
1.9 0.571941 -0.0231087
1.95 0.569791 -0.0430032
}{\numBesselJtable}
\begin{frame}{Example}
Let's take the Bessel equation:\hfill\only<2>{in explicit form}
\begin{overlayarea}{\textwidth}{.2\textheight}
\begin{equation*}
\only<1>{
x^2 y'' + x y' + (x^2 - \nu^2) y = 0
}
\only<2->{
y'' = -\frac{1}{x} y' - \left(1 - \frac{\nu^2}{x^2}\right) y
}
\end{equation*}
\end{overlayarea}
\begin{overlayarea}{\textwidth}{.8\textwidth}
\begin{columns}
\begin{column}{.7\textwidth}
\only<2->{
\resizebox{\textwidth}{!}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=left,
ymin=-0.50,ymax=+0.65,
xmin=+1.00,xmax=+2.00
]
% Bessel J
\addplot+[
ultra thick,
mark=none,
color=blue,
] table [
x index=0,
y index=1
] {\BesselJtable};
\addplot+[
ultra thick,
mark=none,
color=red,
] table [
x index=0,
y index=2
] {\BesselJtable};
\addplot+[
ultra thick,
only marks,
mark=x,
color=blue,
] table [
x index=0,
y index=1
] {\numBesselJtable};
\addplot+[
ultra thick,
only marks,
mark=+,
color=red,
] table [
x index=0,
y index=2
] {\numBesselJtable};
% BesselY1 * (0.44051 / -0.781213)
\only<3>{
\addplot+[
ultra thick,
dashed,
mark=none,
color=blue,
] table [
x index=0,
y index=1
] {\BesselYtable};
\addplot+[
ultra thick,
dashed,
mark=none,
color=red,
] table [
x index=0,
y index=2
] {\BesselYtable};
}
\end{axis}
\end{tikzpicture}
}
}
\end{column}
\begin{column}{.3\textwidth}
\only<2>{
\begin{block}{IVP}
\begin{align*}
y(1) & = 0.44051 \\
y'(1) & = 0.325147
\end{align*}
step size: $h = 0.05$\\
\bigskip
Using implicit Euler method!
\end{block}
}
\only<3>{
\begin{block}{IVP (dashed)}
\begin{align*}
y(1) & = 0.44051 \\
y'(1) & = -0.490276
\end{align*}
step size: $h = 0.05$\\
\bigskip
Using implicit Euler method!
\end{block}
}
\end{column}
\end{columns}
\end{overlayarea}
\end{frame}
\begin{frame}{Boundary Value Problem and the Shooting Method}
\begin{itemize}
\item a \alert{boundary value problem} for a second-order ODE
consists of finding the solution
to the ODE on the interval $[a, b]$ with the following condition
\end{itemize}
\begin{align*}
y(a) & = \alpha &
y(b) & = \beta
\end{align*}
\begin{itemize}
\item we can rewrite it as an initial value problem:
\end{itemize}
\begin{align*}
y(a) & = \alpha &
y'(a) & = \gamma
\end{align*}
\begin{itemize}
\item We can the proceed to solve the IVP as a function of $\gamma$, such
$y(b) = \beta$ is fulfilled.
\item This is called the \alert{shooting method}
\item Practically, we define a root-finding problem:
\end{itemize}
\begin{align*}
0 = F(\gamma) \equiv \beta - y(x = b; \alpha, \gamma)
\end{align*}
\end{frame}
\begin{frame}{Example}
From Stoer, Bulirsch:
\begin{equation*}
w''(t) = \frac{3}{2} w^2(t)
\end{equation*}
with boundary conditions
\begin{align*}
w(0) & = 4 & w(1) & = 1
\end{align*}
\begin{columns}
\begin{column}{.6\textwidth}
{$w'(t) = s$}
\only<1>{
\includegraphics[width=\textwidth]{Shooting_method_trajectories}
}
\only<2->{
\includegraphics[width=\textwidth]{Shooting_method_error}
}
\end{column}
\begin{column}{.4\textwidth}
\only<2->{
\alert{two} solutions:
\begin{enumerate}
\item $w'(0) = -8$
\item $w'(0) \simeq -36$.
\end{enumerate}
}
\end{column}
\end{columns}
\end{frame}
\begin{frame}{Finite Difference Method}
\begin{itemize}
\item totally different approach then previous ones!
\item do not integrate the right-hand side of the ODE!
\item instead, discretize the problem (with $k = 1, \dots, K$)
\end{itemize}
\begin{align*}
x_k & = x_0 + k h & y_k &: \text{independent variables}
\end{align*}
\begin{itemize}
\item discretize the derivatives:
\end{itemize}
\begin{align*}
y'(x_k) & = \frac{y_{k + 1} - y_{k - 1}}{2 h} &
y''(x_k) & = \frac{y_{k + 1} + y_{k - 1} + 2 y_k}{h^2} &
& \dots
\end{align*}
\begin{itemize}
\item boundary value problem now boils down to solving the matrix-valued equation
\end{itemize}
\begin{equation*}
A \cdot \vec{y} = \vec{b}
\end{equation*}
\begin{itemize}
\item $A$ only depends on the ODE!
\item $y^T = (y_1, \dots, y_{K})$
\item $b^T$ to be determined, includes boundary problem in $b_1$ and $b_{K}$
\end{itemize}
\end{frame}
\begin{frame}{Finite Difference Method}
Setup:
\begin{itemize}
\item if the ODE does not depend on $x$, then the matrix $A$ can be prepared early
\item in that case, $A$ only depends on the type of differential operator acting on $y$
\item the size $K$ of the problem can easily be adjusted at run time
\item the inverse of $A$ can be determined at compile time
\end{itemize}
$\Rightarrow$ FDMs can be implemented rather efficiently for ``easy'' differential
equations
\vfill
Convergence:
\begin{equation*}
|y(x_k) - y_k| \leq \frac{y^{(4)}(\xi) h^2}{24} (x_k - a) (b - x_k)
\end{equation*}
for $\xi \in [a, b]$.
\end{frame}
\begin{frame}{Finite Difference Method}
\begin{block}{Example \hfill $K=9$, $h=1/(K+1) = 0.1$}
\begin{equation*}
y'' = 4\qquad y(0) = y_0 = -1\qquad y(1) = y_{10} = +1
\end{equation*}
The analytic solution reads $y(x) = 2x^2 - 1$
\begin{overlayarea}{\textwidth}{.7\textheight}
\only<1>{
\begin{align*}
\frac{1}{h^2}\left(\begin{matrix}
-2 & +1 & 0 & \dots & 0 \\
+1 & -2 & +1 & \ddots & \vdots \\
0 & +1 & \ddots & \ddots & 0 \\
\vdots & \ddots & \ddots & \ddots & +1 \\
0 & \dots & 0 & +1 & -2
\end{matrix}\right)
\cdot
\left(\begin{matrix}
y_1\\
y_2\\
\vdots\\
y_{K-1}\\
y_{K}
\end{matrix}\right)
=
\left(\begin{matrix}
4 + \frac{1}{h^2}\\
4\\
\vdots\\
4\\
4 - \frac{1}{h^2}
\end{matrix}\right)
\end{align*}
}
\only<2>{
\resizebox{!}{.7\textheight}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=left,
ymin=-1.00,ymax=+1.00,
xmin=+0.00,xmax=+1.00,
samples=160
]
\addplot+[
black,
ultra thick,
mark=none
]{
2*x*x - 1
};
\addplot+[
red,
thick,
mark=+,
only marks
] coordinates{
(0.1, -0.98)
(0.2, -0.92)
(0.3, -0.82)
(0.4, -0.68)
(0.5, -0.50)
(0.6, -0.28)
(0.7, -0.02)
(0.8, +0.28)
(0.9, +0.62)
};
\end{axis}
\end{tikzpicture}
}
}
\end{overlayarea}
\end{block}
\end{frame}
\backupbegin
\begin{frame}\frametitle{Backup}
\end{frame}
\backupend
\end{document}
Danny van Dyk