% 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 have 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
\vspace{-\medskipamount}
\begin{equation*}
\sum_{k=0}^{K} \omega^{(K)}_k = 1
\end{equation*}
~\\
\vspace{-1.5\bigskipamount}
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-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<2>{
\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}
}
}
\only<3>{
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-1} \int_{x_{m \cdot N}}^{x_{(m+1)\cdot N}} \hspace{-2em}\mathrm{d}x\, f(x)
= \frac{b - a}{M} \sum_{m=1}^{M-1} \sum_{n=0}^{N} \omega_k f(x_{m\cdot N + n})
= \frac{b - a}{N \cdot M} \hspace{-0.8em}\sum_{k=0}^{K = N \cdot M} \hspace{-0.5em}\Omega_k f(x_k)
\end{align*}
with $x_i = a + (b - a) \frac{i}{N \cdot M}$\,.
}
\end{frame}
\begin{frame}{(Closed) Newton-Cotes Weights and error estimates}
\begin{center}
\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}
\end{center}
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!
\item The iterated weights $\Omega_k$ can be obtained from shuffling and interweaving the single-interval weights $\omega_k$, e.g.:
\end{itemize}
\begin{align*}
\Omega^{(1)}_k & = \lbrace \frac12, 1, 1, \dots, 1, \frac12 \rbrace\\
\Omega^{(2)}_k & = \lbrace \frac16, \frac46, \frac26, \frac46, \frac26, \dots, \frac26, \frac16 \rbrace
\end{align*}
\end{frame}
\begin{frame}{Should we always use equidistant grids?}
Both open and closed Newton-Cotes formula use fixed-size equidistant
grids to interpolate the target function piece by piece.\\
\vfill
Should we always do that?\\
\vfill
Can it be benefitial to use prior knowledge of the target function
in our approximation of the interval? How to incorporate that into
the methods that we know?
\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)
};
\only<2->{
\addplot[thick,blue,only marks,mark=x,mark size=5pt] coordinates {
(-3, 81/82)
(-(11/6), 121/125)
(-(2/3), 4/5)
(1/2, 9/13)
(5/3, 25/26)
(17/6, 289/293)
(4, 144/145)
};
}
\only<3->{
\addplot[thick,red,only marks,mark=x,mark size=5pt] coordinates {
(-3, 81/82)
(-(3/2), 81/85)
(0, 0)
(0, 0)
(2, 36/37)
(4, 144/145)
};
}
\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}
{
\setbeamercolor{block title}{bg=blue}
\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}
}
{
\setbeamercolor{block title}{bg=red}
\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 (e.g.: roots, minima, poles!) 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) p(x)
\end{equation*}
Nomenclature:
\begin{itemize}
\item $p(x)$ is a polynomial of degree $N$.
\item $\omega(x)$ is the weight function: $\omega(x) > 0$, integrable and has finitely many roots in $[a, b]$.
\item $p_i(x)$ are elements of the vector space of orthogonal polynomials under $\omega$:
\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{itemize}
We can then approximate
\begin{equation*}
I \approx \sum_{n=1}^N \omega_n f(\xi_n)
\end{equation*}
where the evaluation points $\xi_n$ with $n=1, \dots, N$ are not yet defined.
\end{frame}
\begin{frame}{General idea cont'd}
What happened? Why is this different from Newton-Cotes?
\begin{itemize}
\item The evaluation points $\xi_n$ are not equidistant anymore!
\item Choose $\omega_n$ and $\xi_n$ such that
\begin{align*}
\tilde{p}(x) & = p(\xi_n) {\color{blue} \left[\prod_{m=1,n\neq m}^{N} \frac{x - \xi_n}{\xi_m - \xi_n}\right]}\\
0 & = \sum_{n=1}^N \omega_n p_i(\xi_n) \qquad\forall i < N
\end{align*}
\end{itemize}
We can achieve this through
\begin{equation*}
\omega_n \equiv \int_{a}^{b} \mathrm{d}x\, \omega(x) \, {\color{blue}\left[\prod_{m=1,n\neq m}^{N} \frac{x - \xi_n}{\xi_m - \xi_n}\right]}
\end{equation*}
and choosing $\xi_n$ as the roots of $p_N(x)$.
\begin{itemize}
\item $\tilde{p}(x)$ corresponds to an interpolating polynomial of
degree $N - 1$ that interpolates $p(x)$ in the $N$ roots of a
(orthogonal) polynomial $p_N(x)$.
\item We now approximate $p(x)$ through $\tilde{p}(x)$.
\end{itemize}
\end{frame}
\begin{frame}{Why is this so benefitial}
\begin{block}{Interpolating polynomial}
\begin{equation*}
\tilde{p}(x) = \sum_n p(\xi_n) \left[\prod_{m=1,n\neq m}^{N} \frac{x - \xi_n}{\xi_m - \xi_n}\right]
\end{equation*}
\end{block}
We can measure the distance between the interpolating polynomial and any other polynomial
of degree $N' = 0, \dots, N - 1 < N$:
\begin{equation*}
\langle p_{N'}, p \rangle = 0 = \int_{a}^{b} \omega(x) p_{N'}(x) p(x)
\end{equation*}
by construction.\\
\vfill
Here $p_{N'}(x) p(x)$ is a polynomial of up to degreen $N' + N \leq 2N - 1$. Therefore
no quadrature formula of degree $< 2N$ exists which yields a better approximation to $I$
than the Gaussian quadrature.
\end{frame}
\begin{frame}{Example}
\begin{columns}
\begin{column}[T]{.5\textwidth}
\begin{itemize}
\item Let {\color{blue}$p(x) = 1 + 3 x + x^2$}, with degree $N=2$.
\item Let $\omega(x) = 1$, which impllies usage of the Legendre polynomials.
\item Let $\xi_n = \lbrace -1 / \sqrt{3}, +1 / \sqrt{3}\rbrace$:
the roots of $P_2$.
\item Interpolate $p(x)$ in the roots:
\begin{equation*}
{\color{red}\tilde{p}(x) = \frac{4}{3} - 3 x}
\end{equation*}
\end{itemize}
\end{column}
\begin{column}[T]{.5\textwidth}
\resizebox{\textwidth}{!}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=left,
ymin=-1.8,ymax=+5,
xmin=-1,xmax=+1,
samples=50
]
\addplot+[thick,blue,mark=none,domain=-1:1,name path=A] {
1 + 3 * x + x^2
};
\addplot+[thick,red,mark=none,domain=-1:1,name path=B] {
4/3 + 3 * x
};
\addplot[gray!50] fill between[of=A and B];
\end{axis}
\end{tikzpicture}
}
\end{column}
\end{columns}
\begin{itemize}
\item The integrals over both $\omega(x) p(x)$ and $\omega(x) \tilde{p}$ coincide: $I = 8/3$.
\item The discrete weights are $\omega_1 = \omega_2 = 1$\\
$\Rightarrow$ $I\approx p(\xi_1) + p(\xi_2) = \left[\frac{4}{3} - \sqrt{3}\right] + \left[\frac{4}{3} + \sqrt{3}\right] = \frac{8}{3}$.
\end{itemize}
\end{frame}
\begin{frame}{List of weight functions and associated polynomials}
\begin{center}
\resizebox{\textwidth}{!}{
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{c|c|cc}
\textbf{interval} $\bm{[a, b]}$ & \textbf{weight} $\bm{\omega(x)}$ & \multicolumn{2}{c}{\textbf{name of polynomials}}\\
\hline
$[-1, 1]$ & $1$ & $P_n(x)$ & Legendre polynomials\\
$[-1, 1]$ & $(1 - x^2)^{-1/2}$ & $T_n(x)$ & Tschebyscheff polynomials\\
$[-1, 1]$ & $(1 - x)^\alpha (1 + x)^\beta$ & $P^{\alpha,\beta}_n(x)$ & Jacobi polynomials\\
$[0, \infty]$ & $\exp(-x)$ & $L_n(x)$ & Laguerre polynomials\\
$[-\infty,\infty]$ & $\exp(-x^2)$ & $H_n(x)$ & Hermite polynomials\\
\end{tabular}
\renewcommand{\arraystretch}{1}
}
\end{center}
You can look up this list, as well as the roots $\xi_n$ and the weights $\omega_n$ in:
\begin{center}
\textit{Handbook of mathematical functions}, Abramowith and Stegun (1964).
\end{center}
\end{frame}
\backupbegin
\begin{frame}\frametitle{Backup}
\end{frame}
\backupend
\end{document}
Danny van Dyk