diff --git a/Lectures_my/NumMet/2016/Lecture7/lecture7.tex b/Lectures_my/NumMet/2016/Lecture7/lecture7.tex index 612d9fd..e2b4626 100644 --- a/Lectures_my/NumMet/2016/Lecture7/lecture7.tex +++ b/Lectures_my/NumMet/2016/Lecture7/lecture7.tex @@ -85,22 +85,105 @@ \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}{Linear approximation} +\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!} - Integrate the interpolating polynomial! + \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} - Newton-Coates! +\begin{frame}{Newton-Coates Quadrature} + Isaac Newton and Roger Cotes had this idea before us! - \begin{equation} - \sum_k \omega_k = 1\,. - \end{equation} + \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} @@ -160,10 +243,12 @@ \end{itemize} \vfill The integral can then be approximated as - \begin{equation} - I \approx \sum_{k=0}^{N \cdot M + 1} \omega_k f(x_k) - \end{equation} - with $x_k = a + (b - a) \frac{k}{N\cdot M}$\,. + \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}{!}{ @@ -181,41 +266,6 @@ \addplot[thick,black,domain=0:+1] { 1 / (1 + 25 * x * x) }; - \only<3>{ - \addplot+[smooth,red,name path=A1,domain=-1:-0.6,mark=none] { - 0.348416 + 0.570136 * x + 0.260181 * x^2 - }; - \addplot+[draw=none,mark=none,name path=B1] {0}; - \addplot+[gray] fill between[of=A1 and B1,soft clip={domain=-1:-0.6}]; - } - \only<4>{ - \addplot+[smooth,red,name path=A2,domain=-0.6:-0.2,mark=none] { - 1. + 3. * x + 2.5 * x^2 - }; - \addplot+[draw=none,mark=none,name path=B2] {0}; - \addplot+[gray] fill between[of=A2 and B2,soft clip={domain=-0.6:-0.2}]; - } - \only<5>{ - \addplot+[smooth,red,name path=A3,domain=-0.2:+0.2,mark=none] { - 1. - 12.5 * x^2 - }; - \addplot+[draw=none,mark=none,name path=B3] {0}; - \addplot+[gray] fill between[of=A3 and B3,soft clip={domain=-0.2:+0.2}]; - } - \only<6>{ - \addplot+[smooth,red,name path=A4,domain=+0.2:+0.6,mark=none] { - 1. - 3. * x + 2.5 * x^2 - }; - \addplot+[draw=none,mark=none,name path=B4] {0}; - \addplot+[gray] fill between[of=A4 and B4,soft clip={domain=+0.2:0.6}]; - } - \only<7>{ - \addplot+[smooth,red,name path=A5,domain=+0.6:+1,mark=none] { - 0.348416 - 0.570136 * x + 0.260181 * x^2 - }; - \addplot+[draw=none,mark=none,name path=B5,domain=+0.6:+1] {0}; - \addplot+[gray] fill between[of=A5 and B5,soft clip={domain=0.6:1}]; - } \addplot[thick,red,only marks,mark=*,mark size=3pt] coordinates { (-1.0, 0.0384615) (-0.8, 0.0588235) @@ -229,20 +279,51 @@ ( 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}{} - \begin{equation} - \omega_k = \begin{cases} - \frac{1}{6} & k = 0 \lor k = N \cdot M + 1\\ - \frac{4}{6} & k\,\text{even}\\ - \frac{2}{6} & k\,\text{odd} - \end{cases} - \end{equation} +\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} @@ -285,7 +366,8 @@ I & = 5.57396 \end{align*} \end{block} - \begin{block}{Entire range, 7 data points} + \vspace{-\smallskipamount} + \begin{block}<2->{Entire range, 7 data points} ~\\ \vspace{-2\bigskipamount} \begin{align*} @@ -293,7 +375,8 @@ \delta_7 & = 12.26\% \end{align*} \end{block} - \begin{block}{Two integration, 5 data points} + \vspace{-\smallskipamount} + \begin{block}<3->{Two integration, 5 data points} ~\\ \vspace{-2\bigskipamount} \begin{align*} @@ -303,6 +386,46 @@ \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