Newer
Older
Lecture_repo / Lectures_my / NumMet / Lecture3 / lecture3.tex
@Danny van Dyk Danny van Dyk on 22 Sep 2016 21 KB Update Lecture 3
% 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}