% 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{tikz,pgfplots}
\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 {Introduction to \\Numerical Methods}
\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 Chrząszcz, Danny van Dyk\\\vspace{-0.1em}\small \href{mailto:mchrzasz@cern.ch}{mchrzasz@cern.ch}, \href{mailto:dany.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}
% \footnotesize\textcolor{gray}{With N. Serra, B. Storaci\\Thanks to the theory support from M. Shaposhnikov, D. Gorbunov}\normalsize\\
\vspace{0.5em}
\textcolor{normal text.fg!50!Comment}{Numerical Methods, \\ 26. September, 2016}
\end{center}
\end{frame}
}
\begin{frame}{Does interpolation always work?}
\begin{block}{Weierstrass Approximation Theorem}
Suppose $f$ is a continuous real-valued function defined on the real interval $[a, b]$.
For every $\varepsilon > 0$, there exists a polynomial $P(x)$ such that for all $x$ in $[a, b]$,
we have $R[P] \equiv \max_{a \leq x \leq b} |f (x) - P(x)| < \varepsilon$, [\dots]
\end{block}
Using interpolation, we can construct a polynomial $P_n(x)$ of degree $n$
that approximates $f$ up to an approximation error of $R_n(x)$:
\begin{equation*}
R_n(x) \equiv \max_{a \leq x \leq b} \left|f(x) - P_n(x)\right|\,.
\end{equation*}
\begin{alert}{Conclusion}
Weierstrass's theorem says that there is \emph{at least one} polynomial
for each choice of the residual error $\varepsilon$. It does not say that either
\begin{itemize}
\item $P_n$ is in the set of polynomial for a given $\varepsilon$, or
\item $R_n \to 0$ if $n \to \infty$!
\end{itemize}
\begin{center}{\color{red} Increasing $n$ might be harmful!}\end{center}
\end{alert}
\end{frame}
\begin{frame}{When Interpolation fails}
So far, you have been confronted with examples in which interpolation
works nicely. Let us now discuss an example, which behaves pathologically.\\[\medskipamount]
\begin{overlayarea}{\textwidth}{6cm}
\begin{columns}
\begin{column}[T]{.5\textwidth}
Consider the classical example by Runge:
\begin{equation*}
f(x) = \left[1 + 25 x^2\right]^{-1}
\end{equation*}
Let us plot
\begin{itemize}
\item the true function $f(x)$
\item<2-> {\color{blue}the interpolating polynomial to degree $4$}
\item<3-> {\color{red}the interpolating polynomial to degree $6$}
\end{itemize}
\end{column}
\begin{column}[T]{.5\textwidth}
\begin{tikzpicture}
[
scale=0.75
]
\begin{axis}[%
samples=500,
xmin=-1,xmax=+1,
ymin=-0.5,ymax=+1
]
\addplot+[black,mark=none]
{1.0 / (1.0 + 25.0 * \x^2)};
\only<2->{
\addplot+[blue,domain=-1:+1,y domain=-1:+1,mark=none]
{1 - (3225 * \x^2) / 754 + (1250 * \x^4) / 377};
}
\only<3->{
\addplot+[red,domain=-1:+1,y domain=-1:+1,mark=none]
{1 - (211600 * \x^2)/24089 + (2019375 * \x^4)/96356 - (1265625 * \x^6)/96356};
}
\end{axis}
\end{tikzpicture}
\end{column}
\end{columns}
\begin{center}\only<4->{$R_n \to \infty$ as $n \to \infty$}\end{center}
\end{overlayarea}
\end{frame}
\begin{frame}{Pieceswise linear/cubic/... interpolations}
Piecewise interpolation of the target function $f$ can be achieved
with relatively small degrees of freedom for each of the interpolating polynomials.
Most popular are linear ($n=1$) or cubic ($n=3$) piecewise interpolations
(``splines''), which do not suffer from the problem that was just described.
\end{frame}
\begin{frame}{Interpolation in more than $D=1$ dimensions}
In lecture 2 we discussed various ways to interpolate a
univariate function $f(x)$. A very nice fact for
polynomial interpolations in $D=1$ dimension is that the
interpolating polynomial is \emph{unique}: For two polynomials
$g(x)$ and $h(x)$ of identical degree $n$ one has
\begin{equation*}
f_i = g(x_i) = h(x_i) \quad\Rightarrow\quad g(x) \equiv h(x)\,,
\end{equation*}
when considering exactly $\rho_n(D=1) = n + 1$ interpolation points $x_i$,
for $i = 1,\dots, n+1$.\\
\medskip
We will now investigate if this still holds for $D > 1$ dimensions.
\end{frame}
\begin{frame}{Interpolation in $D=2$ dimensions}
We will use a bivariate polynomial of degree $n$:
\begin{equation*}
P_n(x, y) = \sum_{i, j}^{i + j \leq n} a_{i, j} x^i y^j\,.
\end{equation*}
This polynomial has exactly $\rho_n(D=2) = \binom{n + 2}{n}$ coefficients,
which need to be computed from the interpolation points.\\
If we evaluate exactly $\rho_n(D=2)$ points, the system of linear equations
has exactly one or zero solutions!
\begin{block}{Examples}
\vspace{-\medskipamount}
\begin{align*}
P_1(x, y) & = a_{0,0} + a_{1,0} x + a_{0,1} y\,, & \binom{1 + 2}{1} = 3\,,\\
P_2(x, y) & = a_{0,0} + a_{1,0} x + a_{0,1} y\,, & \binom{2 + 2}{2} = 6\,,\\
& + a_{1,1} x y + a_{2, 0} x^2 + a_{0, 2} y^2
\end{align*}
\end{block}
\end{frame}
\begin{frame}{Arrangement of points in $D=2$}
The question is now: How do we choose the $\rho_n(D=2)$ points that shall be
interpolated?\\
\medskip
In $D=1$, this was easy: all points were on the $x$ axis. In $D=2$ we need to populate
a plane, and this gives us more freedom.\\
\medskip
How to arrange the interpolation points in the plane? Let's consider two examples
for $D = 2$ dimensions for a $n=1$ polynomial:
\begin{overlayarea}{\textwidth}{5cm}
\only<1>{
\begin{block}{Example \#1}
\begin{columns}
\begin{column}[T]{.35\textwidth}
\begin{center}
\begin{tikzpicture}
\begin{axis}[%
width=3cm,height=3cm,
xmin=0,xmax=1,xtick={0,1},xlabel=$x$,
ymin=0,ymax=1,ytick={0,1},ylabel=$y$
]
\addplot+[blue, only marks, mark options={scale=1.75, fill=blue}] coordinates {
(0.0, 0.0)
(0.5, 0.5)
(1.0, 1.0)
};
\end{axis}
\end{tikzpicture}
\end{center}
\end{column}
\begin{column}[T]{.65\textwidth}
The interpolation equation reads:
\begin{equation*}
\left[\begin{matrix}
1 & 0 & 0 \\
1 & \frac12 & \frac12\\
1 & 1 & 1
\end{matrix}\right]
\cdot
\left[\begin{matrix}
a_{0, 0}\\
a_{1, 0}\\
a_{0, 1}
\end{matrix}\right]
=
\left[\begin{matrix}
f(0, 0)\\
f(\frac12, \frac12)\\
f(1, 1)
\end{matrix}\right]
\end{equation*}
The Vandermonde matrix is singular, and no polynomial exists
that interpolates the blue points.
\end{column}
\end{columns}
\end{block}
}
\only<2>{
\begin{block}{Example \#2}
\begin{columns}
\begin{column}[T]{.35\textwidth}
\begin{center}
\begin{tikzpicture}
\begin{axis}[%
width=3cm,height=3cm,
xmin=0,xmax=1,xtick={0,1},xlabel=$x$,
ymin=0,ymax=1,ytick={0,1},ylabel=$y$
]
\addplot+[red, only marks, mark options={scale=1.75, fill=red}] coordinates {
(0, 0)
(0, 1)
(1, 0)
};
\end{axis}
\end{tikzpicture}
\end{center}
\end{column}
\begin{column}[T]{.65\textwidth}
The interpolation equation reads:
\begin{equation*}
\left[\begin{matrix}
1 & 0 & 0 \\
1 & 1 & 0 \\
1 & 0 & 1
\end{matrix}\right]
\cdot
\left[\begin{matrix}
a_{0, 0}\\
a_{1, 0}\\
a_{0, 1}
\end{matrix}\right]
=
\left[\begin{matrix}
f(0, 0)\\
f(\frac12, \frac12)\\
f(1, 1)
\end{matrix}\right]
\end{equation*}
The Vandermonde matrix is regular, and exactly one polynomial exists
that interpolates the red points.
\end{column}
\end{columns}
\end{block}
}
\end{overlayarea}
\end{frame}
\begin{frame}{Divide and conquer $D=2$}
We had seen that increasing the degree $n$ does not necessarily reduce the
approximation error $R_n$, even in $D=1$ dimensions.\\ \medskip
In $D=2$ this problem can become even more serious.\\ \medskip
It is therefore a good idea to evaluate the function $f$ on a regular grid
in the $(x,y)$ plane. Subequently, one can interpolate within each
cell of the grid. This is the $D=2$ analog to interpolation with splines:
\begin{itemize}
\item linear splines $\to$ bilinear splines
\item cubic splines $\to$ bicubic splines
\end{itemize}
\end{frame}
\begin{frame}{Algorithm for Bilinear Interpolation}
\begin{enumerate}
\item Create a rectilinear grid in the $(x,y)$ plane
\item For each rectangle, map the rectangle to the unit square
\item Evaluate the function $f$ on the four
corners of the (mapped) unit square $P_{0,0}$, $P_{0, 1}$, $P_{1, 1}$, $P_{1, 0}$.
\item Make an ansatz: $f(x, y) = a_{0,0} + a_{1,0} x + a_{0,1} y + a_{1,1} x y$.
\item Solve the interpolation equation
\begin{equation*}
\left[\begin{matrix}
1 & 0 & 0 & 0\\
1 & 1 & 0 & 1\\
1 & 0 & 1 & 1\\
1 & 1 & 1 & 1
\end{matrix}\right]
\cdot
\left[\begin{matrix}
a_{0, 0}\\
a_{1, 0}\\
a_{0, 1}\\
a_{1, 1}
\end{matrix}\right]
=
\left[\begin{matrix}
f(P_{0, 0})\\
f(P_{1, 0})\\
f(P_{0, 1})\\
f(P_{1, 1})
\end{matrix}\right]
\end{equation*}
\item Map the unit square back to the $(x,y)$ plane.
\end{enumerate}
\end{frame}
\begin{frame}{Interpolation in $D$ dimensions}
The situation becomes even more complicated if you go to
$D >2$ dimensions:
\begin{description}
\item[$D=3$] You can generalize to with either trilinear or tricubic splines, on
a rectilinear $3D$ grid. The Vandermonde matrix on the unit cube is now
$\binom{n + 3}{n}^2$.
\item[arbitrary $D$] There is no polynomial interpolation for $D$ dimensions.
\end{description}
\end{frame}
\begin{frame}[shrink]{Extrapolation Basics}
In many numerical applications a common class of problems arises:
In the valuation of a function $f(x)$, we are interested in the value
$f(x_0)$. However, at $x_0$ the function $f(x)$ is numerically instable,
or maybe even ill-posed.\\ \medskip
However, in an environment around $x_0$, $x \approx x_0 + h$, we can evaluate $f$.
Usually, one now discusses $f(h) \equiv f(x_0 + h)$, or rather the limit
$\lim_{h \to 0} f(h)$. [Note: in all generality we can map problems with
limits to either $\infty$ or a finite value to limits to $0$.]\\ \medskip
Interpolation, as discussed previously, can not directly help, since $h = 0$
is not part of the domain of data points. Instead, one can take an interpolation
at finite $h > 0$, $f_\text{int}(h)$, and simply approximate $f(h = 0) \approx f_\text{int}$.
This step of using the interpolation of $f$ outside the domain of data points is
called \emph{extrapolation}.\\ \medskip
To extrapolate an arbitrary function might work very well, but also might
fail spectacularly. In this part of the course, we will briefly discuss
examples of both cases, and what mathematical requirements make extrapolations
work.
\end{frame}
\begin{frame}{Working example}
Consider the function $\cos(x)$. Our task at hand is to estimate
$\cos(0)$ through evaluates of $\cos(h_n)$, with $h_n > 0$ but
$h_n \to 0$ with increasing $n$.\\ \medskip
Let $h_k = 2^{-k}$, and construct interpolating polynomials $P_n$
that interpolate $\cos(x)$ in $\lbrace h_1, \dots, h_{n + 1}\rbrace$. \\ \medskip
\begin{overlayarea}{\textwidth}{5cm}
\begin{columns}
\begin{column}[T]{.5\textwidth}
\only<1>{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0.5, ymax=1.05,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick] {
cos(deg(x))
};
\addplot[blue] {
cos(deg(0.5))
} node [above,pos=1] {$n=0$};
\addplot[black, mark=x] coordinates {
(0.5000, 0.877583)
};
\end{axis}
\end{tikzpicture}
}
\only<2>{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0.5, ymax=1.05,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick] {
cos(deg(x))
};
\addplot[blue] {
1.06 - 0.365 * x
} node [above,pos=1] {$n=1$};
\addplot[black, mark=x] coordinates {
(0.5000, 0.877583)
(0.2500, 0.968912)
};
\end{axis}
\end{tikzpicture}
}
\only<3->{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0.5, ymax=1.05,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick] {
cos(deg(x))
};
\addplot[blue] {
0.99996 + 0.00119746 * x - 0.511201 * x^2 + 0.0385918 * x^3
} node [above,pos=1] {$n=3$};
\addplot[black, mark=x] coordinates {
(0.5000, 0.877583)
(0.2500, 0.968912)
(0.1250, 0.992198)
(0.0625, 0.998048)
};
\end{axis}
\end{tikzpicture}
}
\end{column}
\begin{column}[T]{.5\textwidth}
\only<4->{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
title=$P_n(0)$,
xlabel=$n$
]
\addplot+[black, mark=*, only marks] coordinates {
(0, 0.877583)
(1, 1.06024)
(2, 1.00056)
(3, 0.99996)
(4, 1.00000)
(5, 1.00000)
};
\addplot+[red, mark=none, domain=0:6] { 1 };
\end{axis}
\end{tikzpicture}
}
\end{column}
\end{columns}
\end{overlayarea}
\end{frame}
\begin{frame}{Pathological example}
Consider the example:
\begin{equation*}
f(x) = \exp(-x^{-1/2}) / x^4\,,\quad\text{with}\quad \lim_{x\to 0} f(x) = 0
\end{equation*}
\begin{overlayarea}{\textwidth}{5cm}
\begin{columns}
\begin{column}[T]{.5\textwidth}
\only<1>{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0, ymax=120,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick, samples=250] {
exp(-1 / sqrt(x)) / x^3
};
\addplot[blue] {
1.94493
} node [above,pos=1] {$n=0$};
\addplot[black, mark=x] coordinates {
(0.5000, 1.94493)
};
\end{axis}
\end{tikzpicture}
}
\only<2>{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0, ymax=120,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick, samples=250] {
exp(-1 / sqrt(x)) / x^3
};
\addplot[blue] {
15.378 - 26.8661 * x
} node [above,pos=1] {$n=1$};
\addplot[black, mark=x] coordinates {
(0.5000, 1.94493)
(0.2500, 8.66146)
};
\end{axis}
\end{tikzpicture}
}
\only<3->{
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0, ymax=120,
title=$P_n(x)$,
xlabel=$x$
]
\addplot[black, thick, samples=250] {
exp(-1 / sqrt(x)) / x^3
};
\addplot[blue, samples=50] {
153.618 - 1573.05 * x + 5406.39 * x^2 - 5733.96 * x^3
} node [above,pos=1] {$n=3$};
\addplot[black, mark=x] coordinates {
(0.5000, 1.94493)
(0.2500, 8.66146)
(0.1250, 30.2621)
(0.0625, 75.0209)
};
\end{axis}
\end{tikzpicture}
}
\end{column}
\begin{column}[T]{.5\textwidth}
\only<4->{
\begin{tikzpicture}
[
scale=0.75
]
\begin{axis}[%
ymode=log,
title=$\ln |P_n(0) - 1|$,
xlabel=$n$
]
\addplot+[black, mark=*, only marks] coordinates {
(0, 1.94493)
(1, 15.3780)
(2, 64.0244)
(3, 153.618)
(4, 169.579)
(5, 4.61491)
(6, 94.4333)
(7, 0.598657)
(8, 9.77146)
(9, 1.16233)
};
\addplot+[red, mark=none, domain=0:9] { 1 };
\end{axis}
\end{tikzpicture}
show log of $f(x = 2^{-k})$ for $k=1$ to $10$ (with arbitrary numerical precision!)
}
\end{column}
\end{columns}
\end{overlayarea}
\end{frame}
\begin{frame}{Foundations: Necessary prerequisite for extrapolation}
In order for the extrapolation to work, the function $f(x)$ needs
to fulfill a necessary prerequisite:
\begin{itemize}
\item existence of an asymptotic expansion
\end{itemize}
\begin{block}{Asymptotic Expansion}
A function $f(x)$ can be asymptotically expanded in a sequence $\phi_n(x)$ around $x \simeq x_0$,
\begin{equation*}
f(x) = \sum_n^N a_n \phi_n(x)\,,
\end{equation*}
if
\begin{align*}
f(x) - \sum_n^{N-1} a_n \phi_n(x) = o(\phi_{N-1}(x))\,.
\end{align*}
Example: The Taylor expansion is an asymptotic expansion with $\phi_n = x^n$.
\end{block}
\end{frame}
\begin{frame}{Practical Issues}
The pathological example you saw earlier \emph{does have} an asymptotic expansion!
However, for $n \leq 10$ the extrapolation did not yet show convergence.\\ \medskip
\begin{columns}
\begin{column}[T]{.5\textwidth}
The problem here is the starting point, and the existence of a bump between the starting
point $h_1 = 1/2$ and the point of interest ($h=0$).
\end{column}
\begin{column}[T]{.5\textwidth}
\begin{tikzpicture}
[
scale=0.7
]
\begin{axis}[%
domain=0:1,
ymin=0, ymax=120,
title=$f(x)$,
xlabel=$x$
]
\addplot[black, thick, samples=250] {
exp(-1 / sqrt(x)) / x^3
};
\end{axis}
\end{tikzpicture}
\end{column}
\end{columns}
\end{frame}
\backupbegin
\begin{frame}\frametitle{Backup}
\end{frame}
\backupend
\end{document}
Danny van Dyk