% 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}
\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 {Numerical Integration}
		\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, \\ 26. September, 2016}
\end{center}
\end{frame}
}

\begin{frame}{Plan for today}
    \begin{itemize}
        \item \alert{General problem of estimating the integral of a continuous function $f(x)$ on a finite support}\\
        \vfill
        \item \alert{Specific problems in which properties of the integrand can be used to our advantage}\\
            What properties of $f(x)$ can we use make our live easier?
    \end{itemize}
\end{frame}

\begin{frame}{Setup}
    Commonly, the integral $I$ of a function $f(x)$ over finite support $a \leq x \leq b$
    is introduced via the area below the curve $\gamma: (x, f(x))$. This isually called
    a \alert{Riemann}-type integral.\\
    \vfill
    For this lecture we will stick to this type of integrals. However,
    note that there are other measures that can be used to determine an integral (see \alert{Lebesque}-type
    integrals).\\
    \vfill
    For this first lecture, we will stick to one-dimensional, definite integrals:
    \begin{equation*}
        I \equiv \int_a^b \mathrm{d}x\, f(x)
    \end{equation*}
    Moreover, we will focus on integrals of smooth, continuous functions.
\end{frame}

\begin{frame}{Constant approximation}
    Usually, the first step in introducing students to the concept of integration
    is the constant approximation of the integral through summation of $K$ rectangles
    of equal widths:\\
    \begin{equation*}
        I^\text{rect}_K \equiv \frac{(b - a)}{K} \sum_{k=1}^K f(x_k)\,\qquad x_k = a + (b - a) \frac{k}{K}
    \end{equation*}
    Variations of this use different heights of the rectangles (left-hand, right-hand, central).\\
    \vfill
    In the limit of infinitely many rectangles $K$ or infinitely small steps $(b - a) / K$
    the approximation can be expected to converge toward the true integral $I$:
    \begin{equation*}
        \lim_{K\to \infty} I^\text{rect}_K = I
    \end{equation*}
\end{frame}

\begin{frame}{More sophisticated approximations}
    We can now try to be more sophisticated. Let's approcximate the function not by
    a sequence of rectangles, but by a sequence of trapezoids!
    \begin{equation*}
        I^\text{trap}_K
            \equiv \frac{(b - a)}{K} \sum_{k=0}^{K - 1} \frac{f(x_k) + f(x_{k+1})}{2}\,,\qquad x_k = a + (b - a)\frac{k}{K}
    \end{equation*}

    Compare the two approximations piece by piece:
    \begin{align*}
        \frac{K\, I^\text{rect}_K}{b - a}
            & = 0 &
            & + f(x_1) &
            & + f(x_2) &
            & + \dots &
            & + f(x_{K-1}) &
            & + f(x_K)\\
        \frac{K\, I^\text{trap}_K}{b - a}
            & = \frac12 f(x_0) &
            & + f(x_1) &
            & + f(x_2) &
            & + \dots &
            & + f(x_{K-1}) &
            & + \frac12 f(x_K)\\
    \end{align*}
    \vfill
    This looks familiar!
\end{frame}

\begin{frame}{Wait a minute, I know this!}
    \begin{itemize}
        \item For the rectangles, we use a constant approximation of $f$ within a small interval.
        \item For the trapezoid, we use a linear approimation of $f$ within a small interval.
    \end{itemize}
    We covered this: Let's interpolate the function $f$ on an equidistant grid with grid size $(b - a) / K$,
    and then integrate the interpolating polynomial!
    \vfill
    Polynomial interpolation yields one unique function, regardless of its representation. For
    this part, the representation via Lagrange basis polynomials $\ell_k(x)$ is benefitial:
    \begin{equation*}
        L(x) \equiv \sum_k f(x_k) \ell_k(x)
    \end{equation*}
    We can then proceed to integrate:
    \begin{equation*}
        I \approx \int_a^b \mathrm{d}x\, L(x) = \sum_k f(x_k) \int_a^b \mathrm{d}x\,\ell_k(x)
    \end{equation*}
\end{frame}

\begin{frame}{Newton-Coates Quadrature}
    Isaac Newton and Roger Cotes had this idea before us!

    \begin{block}{\alert{Closed} Newton-Cotes formulas}
    \begin{equation*}
        I^\text{NC}_{K} = \sum_{k=0}^{K} \omega^{(K)}_k f(x_k)\,,\qquad x_k = a + (b - a) \frac{k}{K}
    \end{equation*}
    where the \alert{weights} $\omega^{(K)}_k$ depend on the index $k$ and the degree of the interpolating
    polynomial $K$.
    \end{block}
    The \alert{open} Newton-Cotes formulas are obtained by interpolating equidistantly only \alert{within}
    the support, and not on its boundaries.
    \vfill
    \begin{block}{Consistency check}
    Any given (closed) Newton-Coates formula must fulfill
    \begin{equation*}
        \sum_{k=0}^{K} \omega^{(K)}_k = 1
    \end{equation*}
    To be proven in the exercises.
    \end{block}
\end{frame}

\begin{frame}{Pathological example 1}
    \resizebox{.9\textwidth}{!}{
    \begin{tikzpicture}
        \begin{axis}[%
            axis x line=center,
            axis y line=left,
            ymin=0.0,ymax=1.05,
            xmin=-1,xmax=+1,
            samples=160
        ]
            \addplot[thick,black,domain=-1:0] {
                1 / (1 + 25 * x * x)
            };
            \addplot[thick,black,domain=0:+1] {
                1 / (1 + 25 * x * x)
            };
            \addplot[thick,red,  domain=-1:0] {
                1.0 - 16.8552 * x^2 + 123.36 * x^4 - 381.434 * x^6 + 494.91 * x^8 - 220.942 * x^(10)
            };
            \addplot[thick,red,  domain=0:+1] {
                1.0 - 16.8552 * x^2 + 123.36 * x^4 - 381.434 * x^6 + 494.91 * x^8 - 220.942 * x^(10)
            };
            \addplot[thick,red,only marks,mark=*,mark size=3pt] coordinates {
                (-1.0, 0.0384615)
                (-0.8, 0.0588235)
                (-0.6, 0.1)
                (-0.4, 0.2)
                (-0.2, 0.5)
                ( 0.0, 1.0)
                ( 0.2, 0.5)
                ( 0.4, 0.2)
                ( 0.6, 0.1)
                ( 0.8, 0.0588235)
                ( 1.0, 0.0384615)
            };
        \end{axis}
    \end{tikzpicture}
    }
\end{frame}

\begin{frame}{Conclusion 1: We know how to fix it!}
    \only<1>{
    When interpolating a function, we saw that using splines fixes
    the problem of oscillating interpolating polynomials.\\
    \vfill

    We can now use this to stabilise our interpolation formula as well.\\
    \vfill

    \begin{itemize}
        \item Let $M$ be the number of intervals in which we want to integrate.
        \item Let $N$ be the degree of the interpolating polynomial in each interval.
        \item Enforce that the interpolating function is continous, but not differentiable
            at the interval boundaries.
    \end{itemize}
    \vfill
    The integral can then be approximated as
    \begin{align*}
        I = \sum_{m=1}^{M} \int_{x_{m \cdot N}}^{x_{(m+1)\cdot N}} \mathrm{d}x\, f(x)
          = \sum_{m=1}^{M} \sum_{n=0}^{N} \omega_k f(x_{m\cdot N + n})
          = \sum_{k=0}^{K = N \cdot M} \Omega_k f(x_k)
    \end{align*}
    with $x_k = a + (b - a) \frac{k}{K\cdot M}$\,.
    }
    \only<2->{
    \resizebox{.9\textwidth}{!}{
    \begin{tikzpicture}
        \begin{axis}[%
            axis x line=center,
            axis y line=left,
            ymin=0.0,ymax=1.05,
            xmin=-1,xmax=+1,
            samples=160
        ]
            \addplot[thick,black,domain=-1:0] {
                1 / (1 + 25 * x * x)
            };
            \addplot[thick,black,domain=0:+1] {
                1 / (1 + 25 * x * x)
            };
            \addplot[thick,red,only marks,mark=*,mark size=3pt] coordinates {
                (-1.0, 0.0384615)
                (-0.8, 0.0588235)
                (-0.6, 0.1)
                (-0.4, 0.2)
                (-0.2, 0.5)
                ( 0.0, 1.0)
                ( 0.2, 0.5)
                ( 0.4, 0.2)
                ( 0.6, 0.1)
                ( 0.8, 0.0588235)
                ( 1.0, 0.0384615)
            };
            \only<3->{
                \addplot+[smooth,red,solid,name path=A1,domain=-1:-0.6,mark=none] {
                    0.348416 + 0.570136 * x + 0.260181 * x^2
                };
                \addplot+[smooth,red,solid,name path=A2,domain=-0.6:-0.2,mark=none] {
                    1. + 3. * x + 2.5 * x^2
                };
                \addplot+[smooth,red,solid,name path=A3,domain=-0.2:+0.2,mark=none] {
                    1. - 12.5 * x^2
                };
                \addplot+[smooth,red,solid,name path=A4,domain=+0.2:+0.6,mark=none] {
                    1. - 3. * x + 2.5 * x^2
                };
                \addplot+[smooth,red,solid,name path=A5,domain=+0.6:+1,mark=none] {
                    0.348416 - 0.570136 * x + 0.260181 * x^2
                };
            }
        \end{axis}
    \end{tikzpicture}
    }
    }
\end{frame}

\begin{frame}{(Closed) Newton-Cotes Weights and error estimates}
    %

    {
    \centering
    \renewcommand{\arraystretch}{1.2}
    \begin{tabular}{c|c|c}
        $\bm{N}$ & $\bm{\omega^{(N)}_k}$                             & \textbf{approx. error}\\
        \hline
        $1$      & $\frac{1}{2}\, \frac{1}{2}$                       & $\sim (b - a) f''(\xi)$\\
        $2$      & $\frac{1}{6}\, \frac{4}{6}\, \frac{1}{6}$         & $\sim (b - a) f^{(4)}(\xi)$\\
        $3$      & $\frac18\, \frac38\, \frac38\, \frac18$           & $\sim (b - a) f^{(4)}(\xi)$\\
        $\vdots$ & $\vdots$                                          & $\vdots$\\
    \end{tabular}
    \renewcommand{\arraystretch}{1}
    }

    Comments and observations:
    \begin{itemize}
        \item $\xi \in [a, b]$
        \item While numerical factors change, the approximation error scales only with an even derivative of the integrand!
    \end{itemize}
\end{frame}

\newcommand{\eps}{\varepsilon}
\begin{frame}{Pathological example 2}
    \begin{equation*}
        I \equiv \int_{-3}^{+4} \mathrm{d}x\, \left[1 - \frac{\eps}{x^2 + \eps}\right]
        \, \quad \text{with }\eps = \frac{1}{9}
    \end{equation*}
    \begin{columns}
    \begin{column}{.5\textwidth}
    \resizebox{\textwidth}{!}{
    \begin{tikzpicture}
        \begin{axis}[%
            axis x line=center,
            axis y line=left,
            ymin=0.0,ymax=1.05,
            xmin=-3,xmax=+4,
            samples=160
        ]
            \addplot[thick,black,domain=-3:0] {
                1 - 1/9 / (x^2 + 1/9)
            };
            \addplot[thick,black,domain=0:+4] {
                1 - 1/9 / (x^2 + 1/9)
            };
            \addplot[thick,blue,only marks,mark=+,mark size=5pt] coordinates {
                ( 1, 0.54030 )
            };
        \end{axis}
    \end{tikzpicture}
    }
    \begin{equation*}
    \end{equation*}
    \end{column}
    \begin{column}{.5\textwidth}
        \begin{block}{Exact result}
            ~\\
            \vspace{-2\bigskipamount}
            \begin{align*}
                I      & = 5.57396
            \end{align*}
        \end{block}
        \vspace{-\smallskipamount}
        \begin{block}<2->{Entire range, 7 data points}
            ~\\
            \vspace{-2\bigskipamount}
            \begin{align*}
                I_7      & = 6.25746 &
                \delta_7 & = 12.26\%
            \end{align*}
        \end{block}
        \vspace{-\smallskipamount}
        \begin{block}<3->{Two integration, 5 data points}
            ~\\
            \vspace{-2\bigskipamount}
            \begin{align*}
                I_5      & = 5.65645 &
                \delta_5 & = 1.47\%
            \end{align*}
        \end{block}
    \end{column}
    \end{columns}
    \only<4->{
    It is benefitial to arrange the integration intervals around known special points of the integrand!
    }
\end{frame}

\begin{frame}{Rationale for Gauss quadrature}
    Given that we ended up (again) using polynomial interpolation to numerically solve a problem,
    we can ask ourselves the following question:\\
    \vfill
    \alert{Is there a method to numerically obtain the exact integral of a polynomial of degree $2n - 1$?}
    \vfill
    Carl Friedrich Gau\ss{} answered in the affirmative!
\end{frame}

\begin{frame}{General idea}
    We restrict ourselves to the interval $[a, b]$, and to integrals of the type
    \begin{equation*}
        I = \int_{a}^{b} \mathrm{d}x\, \omega(x) f(x)
    \end{equation*}
    Nomenclature:
    \begin{itemize}
        \item $\omega(x)$ is the weight function: $\omega(x) > 0$, integrable and has finitely many roots.
        \item $p_N(x)$ is a polynomial of degree $N$ over $x$.
        \item $\xi_{n}$, with $n=1,\dots,N$ are the roots of $p_N(x)$.
    \end{itemize}
    Orthogonality relations apply:
    \begin{equation*}
        \langle p_i(x), p_j(x)\rangle \equiv \int_{a}^{b} \mathrm{d}x\, \omega(x) p_i(x) p_j(x) = \delta_{ij}
    \end{equation*}
\end{frame}

\begin{frame}{General idea cont'd}
    We can then approximate
    \begin{equation*}
        I \approx \sum_{n=1}^N \alpha_n f(\xi_n)
    \end{equation*}
    where the coefficients $\alpha_n$ can be determined as
    \begin{equation*}
        \alpha_n = \int_{a}^{b} \mathrm{d}x\, \omega(x) \prod_{m=1,n\neq m}^{N} \frac{x - \xi_n}{\xi_m - \xi_n}
    \end{equation*}
\end{frame}

\backupbegin

\begin{frame}\frametitle{Backup}

\end{frame}

\backupend
\end{document}