% 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 {Interpolation in $D>2$\\ and Extrapolation}
		\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: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}{Plan for today}
    \begin{itemize}
        \item \alert{When does interpolation fail?}\\
            Marcin showed you how to interpolate in $D=1$. What are the pitfalls? And can we circumvent them?
        \vfill
        \item \alert{How to interpolate a function of $D>1$ variables?}\\
            Is there are sensible generalization of interpolation to multivariate functions? If yes,
            let's apply what we learned before to $D=2$.
        \vfill
        \item \alert{What happens when I use my interpolation beyond the domain of data?}\\
            Why would I do that? When can I do that? How well does it work, and how can I improve it?
    \end{itemize}
\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$:
    \begin{equation*}
        R_n \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 \emph{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.\\

\begin{overlayarea}{\textwidth}{7cm}
\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}
\only<2->{for equidistant interpolation points on $[-1, +1]$.}
\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}
    \begin{center}\only<4->{$R_n \to \infty$ as $n \to \infty$}\end{center}
\end{column}
\end{columns}
\end{overlayarea}
\end{frame}

\begin{frame}{Pieceswise linear/cubic/... interpolations}
    \begin{columns}
    \begin{column}{.6\textwidth}
    This pathological behaviour can be avoided by piecewise interpolation with a small value of
    the polynomial degree $n$.
    Most popular are linear ($n=1$) or cubic ($n=3$) piecewise interpolations
    (``splines''), which do not suffer from the problem that was just described.\\ \medskip
    \end{column}
    \begin{column}{.4\textwidth}
        \includegraphics[width=.9\textwidth]{runge.pdf}
    \end{column}
    \end{columns}
    \medskip
    The concrete approaches have been discussed last lecture by Marcin.
\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}[shrink]{Algorithm for Bilinear Interpolation}
    \begin{columns}
    \begin{column}[T]{.8\textwidth}
    \begin{enumerate}
        \item<1-> Create a  rectilinear grid in the $(x,y)$ plane
        \item<2-> For each rectangle, map the rectangle to the unit square
        \item<3-> Evaluate the function $f$ on the four
            corners of the (mapped) unit square ${\color{purple}Q_{0,0}}$,
            ${\color{red}Q_{1, 0}}$, ${\color{orange}Q_{1, 1}}$, ${\color{Gold}Q_{0, 1}}$.
        \item<4-> Make an ansatz: $P(x, y) = a_{0,0} + a_{1,0} x + a_{0,1} y + a_{1,1} x y$.
    \end{enumerate}
    \end{column}
    \begin{column}[T]{.2\textwidth}
        \begin{overlayarea}{\textwidth}{.3\textheight}
        \resizebox{\textwidth}{!}{
        \begin{tikzpicture}
            \draw[black,thick] (0,0) -- (2,0) -- (2,1) -- (0,1) -- (0,0);
            \draw[black,thick] (2,0) -- (4,0) -- (4,1) -- (2,1) -- (2,0);
            \draw[black,thick] (0,1) -- (2,1) -- (2,2) -- (0,2) -- (0,1);
            \draw[black,thick] (2,1) -- (4,1) -- (4,2) -- (2,2) -- (2,1);
            \only<2->{
                \fill[gray,thick] (0,1) -- (2,1) -- (2,2) -- (0,2) -- (0,1);
            }
            \only<3->{
                \fill[purple] (0,1) circle [radius=2pt];
                \fill[red]    (2,1) circle [radius=2pt];
                \fill[orange] (2,2) circle [radius=2pt];
                \fill[Gold]   (0,2) circle [radius=2pt];
            }
        \end{tikzpicture}
        }\\
        \only<2->{
        \resizebox{\textwidth}{!}{
        \begin{tikzpicture}
            \draw[black,thick] (0,0) -- (1,0) -- (1,1) -- (0,1) -- (0,0);
            \only<3->{
                \fill[purple] (0,0) circle [radius=2pt];
                \fill[red]    (1,0) circle [radius=2pt];
                \fill[orange] (1,1) circle [radius=2pt];
                \fill[Gold]   (0,1) circle [radius=2pt];
            }
        \end{tikzpicture}
        }
        }
        \end{overlayarea}
    \end{column}
    \end{columns}
    \begin{enumerate}
        \setcounter{enumi}{4}
        \item<5-> 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(Q_{0, 0})\\
                f(Q_{1, 0})\\
                f(Q_{0, 1})\\
                f(Q_{1, 1})
            \end{matrix}\right]
        \end{equation*}
            We only need to invert the Vandermonde matrix once!
        \item<6-> Map the result back from the unit square back to the grid piece in $(x,y)$ plane.
        \item<7-> Jump back to \#2.
    \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:
    Concerning a function $f(x)$ that is expensive to evaluate, we are interested in the value
    $f(x_0)$. However, at $x_0$ the function $f(x)$ is numerically instable,
    or the problem maybe even ill-posed.\\ \medskip

    Sometime there is an environment around $x_0$, $x \approx x_0 + h$, in which 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 $h \to \infty$ or $h \to \text{finite value}$ onto limits $h \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}(h = 0)$.
    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 evaluations of $\cos(h_k)$, with $h_k > 0$ but
$h_k \to 0$ with increasing $k$.\\ \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*}
We use the same procedure as before:

\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)
            };
        \end{axis}
    \end{tikzpicture}
    }
    \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 are
        \begin{itemize}
            \item the choice of the starting point,
            \item the existence of a bump between the starting
                  point $h_1 = 1/2$ and the point of interest ($h=0$).
        \end{itemize}
    \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}

\begin{frame}{Speeding up extrapolations}
    Let's assume you have carried out an extrapolation to some point $x_0$
    for several different values of the degree $n$:
    \begin{align*}
        P_n(x) & :\,\text{interpolating polynomial} &
        p_n    & \equiv P_n(x_0)\,.
    \end{align*}
    The $p_n$ are a sequence that (hopefully) converges to the value you seek.\\

    For definiteness, let's use the previous example of $f(x) = \cos(x)$. The first
    few elements of the sequence $p_n$ were:
    \begin{align*}
        p_0 & = 0.877583 &
        p_1 & = 1.06024  &
        p_2 & = 1.00056  \\
        p_3 & = 0.99996  &
        p_4 & = 1.00000
    \end{align*}

    Assuming that any sequence $p_n$ is convergent, it would be nice to have a way
    to accelerate the sequence without any additional (potentially costly!) evaluations
    of the function $f(x)$!
\end{frame}

\begin{frame}{Delta-Square rule by Aitkens}
    Define a new sequence $q_n$ as follows:
    \begin{equation*}
        q_n \equiv p_n - \frac{\left[\Delta(p)_{n}\right]^2}{\Delta^2(p)_{n}}\,,
    \end{equation*}
    where
    \begin{align*}
        \Delta(p)_k   & = p_{k + 1} - p_k\,, &
        \Delta^2(p)_k & = \Delta(p)_{k + 1} - \Delta(p)_k\,.
    \end{align*}
    This could be written simpler, but in order to avoid numerical instabilities (floating point numbers!)
    the version written here is preferable.\\ \vfill

    [This method has further applications which we will revisit next lecture.]
\end{frame}

\begin{frame}{Accelerating our example}
    \begin{columns}
    \begin{column}{.5\textwidth}
    \begin{center}
        \begin{tikzpicture}
            \node[anchor=south] (0, 0) {
            \begin{tabular}{lrr}
                $n$  & $p_n$        & $q_{n - 2}$\\
                0    & 0.87758      & ---      \\
                1    & 1.06024      & ---      \\
                2    & 1.00056      & 0.999954 \\
                3    & 0.99996      & 0.999998 \\
                4    & 1.00000      & \dots    \\
            \end{tabular}
            };
            \only<2>{
                \draw[blue,thick] (-0.40, +1.20) rectangle (+0.05, +1.55);
                \draw[blue,thick] (+1.30, +1.20) rectangle (+1.90, +1.55);
                \draw[draw=none, fill=blue, fill opacity=0.25] (-1.30, +1.20) rectangle (+0.05, 2.55);
            }
            \only<3>{
                \draw[red, thick] (-0.18, +0.70) rectangle (+0.05, +1.05);
                \draw[red, thick] (+1.50, +0.70) rectangle (+1.90, +1.05);
                \draw[draw=none, fill=red,  fill opacity=0.25] (-1.30, +0.70) rectangle (+0.05, 2.05);
            }
        \end{tikzpicture}
    \end{center}
    \end{column}
    \begin{column}{.5\textwidth}
    \begin{itemize}
        \item compare $p_n$ with $q_{n - 2}$ to keep things fair\\
            i.e.: $p_2$ vs $q_0$, $p_3$ vs $q_1$ etc.
        \item<2-> {\color{blue}$q_0(p_0, p_1, p_2)$ is competitive with $p_3$!}
        \item<3-> {\color{red} $q_1(p_1, p_2, p_3)$ is competitive with $p_4$!}
    \end{itemize}
    \end{column}
    \end{columns}
\end{frame}

\backupbegin   

\begin{frame}\frametitle{Backup}


\end{frame}




\backupend			
\end{document}
