diff --git a/Lectures_my/NumMet/2016/Lecture7/lecture7.tex b/Lectures_my/NumMet/2016/Lecture7/lecture7.tex index e2b4626..ef28687 100644 --- a/Lectures_my/NumMet/2016/Lecture7/lecture7.tex +++ b/Lectures_my/NumMet/2016/Lecture7/lecture7.tex @@ -150,7 +150,7 @@ \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$, + 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 @@ -179,9 +179,12 @@ \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} @@ -227,30 +230,7 @@ \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->{ + \only<1-2>{ \resizebox{.9\textwidth}{!}{ \begin{tikzpicture} \begin{axis}[% @@ -279,7 +259,7 @@ ( 0.8, 0.0588235) ( 1.0, 0.0384615) }; - \only<3->{ + \only<2>{ \addplot+[smooth,red,solid,name path=A1,domain=-1:-0.6,mark=none] { 0.348416 + 0.570136 * x + 0.260181 * x^2 }; @@ -300,13 +280,33 @@ \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} - % - - { - \centering + \begin{center} \renewcommand{\arraystretch}{1.2} \begin{tabular}{c|c|c} $\bm{N}$ & $\bm{\omega^{(N)}_k}$ & \textbf{approx. error}\\ @@ -317,13 +317,29 @@ $\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} @@ -349,9 +365,27 @@ \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 ) + \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} } @@ -367,6 +401,8 @@ \end{align*} \end{block} \vspace{-\smallskipamount} + { + \setbeamercolor{block title}{bg=blue} \begin{block}<2->{Entire range, 7 data points} ~\\ \vspace{-2\bigskipamount} @@ -375,6 +411,9 @@ \delta_7 & = 12.26\% \end{align*} \end{block} + } + { + \setbeamercolor{block title}{bg=red} \vspace{-\smallskipamount} \begin{block}<3->{Two integration, 5 data points} ~\\ @@ -384,10 +423,11 @@ \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! + It is benefitial to arrange the integration intervals around known special points (e.g.: roots, minima, poles!) of the integrand! } \end{frame} @@ -395,7 +435,7 @@ 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$?} + \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} @@ -403,29 +443,128 @@ \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) + I = \int_{a}^{b} \mathrm{d}x\, \omega(x) p(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)$. + \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} - Orthogonality relations apply: + We can then approximate \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} + 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} - We can then approximate + 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*} - I \approx \sum_{n=1}^N \alpha_n f(\xi_n) + \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*} - where the coefficients $\alpha_n$ can be determined as + 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*} - \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} + \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 diff --git a/Lectures_my/NumMet/2016/Lecture7/subdivision.nb b/Lectures_my/NumMet/2016/Lecture7/subdivision.nb index 1388964..bdf7dc1 100644 --- a/Lectures_my/NumMet/2016/Lecture7/subdivision.nb +++ b/Lectures_my/NumMet/2016/Lecture7/subdivision.nb @@ -10,10 +10,10 @@ NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 158, 7] -NotebookDataLength[ 22366, 528] -NotebookOptionsPosition[ 21255, 486] -NotebookOutlinePosition[ 21592, 501] -CellTagsIndexPosition[ 21549, 498] +NotebookDataLength[ 22574, 528] +NotebookOptionsPosition[ 21461, 486] +NotebookOutlinePosition[ 21798, 501] +CellTagsIndexPosition[ 21755, 498] WindowFrame->Normal*) (* Beginning of Notebook Content *) @@ -326,7 +326,7 @@ CellChangeTimes->{{3.685522731198388*^9, 3.6855228125078278`*^9}, { 3.685522862545142*^9, 3.685522874429077*^9}, {3.685523073648672*^9, 3.685523196670637*^9}, 3.685523697644681*^9, 3.685523881702211*^9, { - 3.685523920312171*^9, 3.6855239436288424`*^9}}] + 3.685523920312171*^9, 3.6855239436288424`*^9}, 3.685684006861432*^9}] }, Open ]], Cell[CellGroupData[{ @@ -353,11 +353,11 @@ FractionBox[ RowBox[{"ArcTan", "[", "8", "]"}], "2"]}]], "Output", CellChangeTimes->{{3.685523511725808*^9, 3.685523522484685*^9}, - 3.68552389016968*^9, 3.685523946035872*^9}], + 3.68552389016968*^9, 3.685523946035872*^9, 3.6856840072559958`*^9}], Cell[BoxData["5.573955509185797`"], "Output", CellChangeTimes->{{3.685523511725808*^9, 3.685523522484685*^9}, - 3.68552389016968*^9, 3.68552394603734*^9}] + 3.68552389016968*^9, 3.685523946035872*^9, 3.685684007257372*^9}] }, Open ]], Cell[CellGroupData[{ @@ -397,7 +397,7 @@ Cell[BoxData["6.2574588650388545`"], "Output", CellChangeTimes->{{3.685523481304573*^9, 3.685523485045783*^9}, 3.685523531318762*^9, {3.685523804103651*^9, 3.685523819289815*^9}, - 3.6855238991612787`*^9, 3.685523947746043*^9}] + 3.6855238991612787`*^9, 3.685523947746043*^9, 3.685684007377079*^9}] }, Open ]], Cell[CellGroupData[{ @@ -410,7 +410,7 @@ CellChangeTimes->{{3.685523996844348*^9, 3.685523998276311*^9}}], Cell[BoxData["0.12262447282305239`"], "Output", - CellChangeTimes->{3.685523998952448*^9}] + CellChangeTimes->{3.685523998952448*^9, 3.68568400774715*^9}] }, Open ]], Cell[CellGroupData[{ @@ -458,17 +458,17 @@ Cell[BoxData["2.399784791965567`"], "Output", CellChangeTimes->{ 3.685523772217701*^9, {3.685523833015215*^9, 3.6855238572254143`*^9}, - 3.685523902925892*^9, 3.6855239501498938`*^9}], + 3.685523902925892*^9, 3.6855239501498938`*^9, 3.685684007800902*^9}], Cell[BoxData["3.256663560111836`"], "Output", CellChangeTimes->{ 3.685523772217701*^9, {3.685523833015215*^9, 3.6855238572254143`*^9}, - 3.685523902925892*^9, 3.685523950151155*^9}], + 3.685523902925892*^9, 3.6855239501498938`*^9, 3.685684007804133*^9}], Cell[BoxData["5.656448352077403`"], "Output", CellChangeTimes->{ 3.685523772217701*^9, {3.685523833015215*^9, 3.6855238572254143`*^9}, - 3.685523902925892*^9, 3.685523950152231*^9}] + 3.685523902925892*^9, 3.6855239501498938`*^9, 3.685684007804735*^9}] }, Open ]], Cell[CellGroupData[{ @@ -481,7 +481,7 @@ CellChangeTimes->{{3.68552397854193*^9, 3.685523989548272*^9}}], Cell[BoxData["0.014799695253336416`"], "Output", - CellChangeTimes->{3.685523989813211*^9}] + CellChangeTimes->{3.685523989813211*^9, 3.6856840079160843`*^9}] }, Open ]] }, WindowSize->{1054, 1179}, @@ -505,30 +505,30 @@ Cell[1879, 56, 1181, 32, 187, "Input"], Cell[CellGroupData[{ Cell[3085, 92, 1343, 29, 209, "Input"], -Cell[4431, 123, 11863, 205, 241, "Output"] +Cell[4431, 123, 11885, 205, 241, "Output"] }, Open ]], Cell[CellGroupData[{ -Cell[16331, 333, 503, 13, 55, "Input"], -Cell[16837, 348, 283, 7, 49, "Output"], -Cell[17123, 357, 157, 2, 32, "Output"] +Cell[16353, 333, 503, 13, 55, "Input"], +Cell[16859, 348, 307, 7, 49, "Output"], +Cell[17169, 357, 180, 2, 32, "Output"] }, Open ]], Cell[CellGroupData[{ -Cell[17317, 364, 1041, 30, 55, "Input"], -Cell[18361, 396, 234, 3, 32, "Output"] +Cell[17386, 364, 1041, 30, 55, "Input"], +Cell[18430, 396, 256, 3, 32, "Output"] }, Open ]], Cell[CellGroupData[{ -Cell[18632, 404, 197, 5, 32, "Input"], -Cell[18832, 411, 89, 1, 32, "Output"] +Cell[18723, 404, 197, 5, 32, "Input"], +Cell[18923, 411, 110, 1, 32, "Output"] }, Open ]], Cell[CellGroupData[{ -Cell[18958, 417, 1387, 38, 99, "Input"], -Cell[20348, 457, 188, 3, 32, "Output"], -Cell[20539, 462, 186, 3, 32, "Output"], -Cell[20728, 467, 186, 3, 32, "Output"] +Cell[19070, 417, 1387, 38, 99, "Input"], +Cell[20460, 457, 210, 3, 32, "Output"], +Cell[20673, 462, 210, 3, 32, "Output"], +Cell[20886, 467, 210, 3, 32, "Output"] }, Open ]], Cell[CellGroupData[{ -Cell[20951, 475, 195, 5, 32, "Input"], -Cell[21149, 482, 90, 1, 32, "Output"] +Cell[21133, 475, 195, 5, 32, "Input"], +Cell[21331, 482, 114, 1, 32, "Output"] }, Open ]] } ]