eth-summaries/semester3/numcs/parts/02_quadrature/00_introduction.tex

% ┌                                                ┐
% │     AUTHOR: Janis Hutz<info@janishutz.com>     │
% └                                                ┘

\setcounter{subsection}{2}
\subsection{Grundbegriffe und -Ideen}
Es ist oft nicht möglich oder sinnvoll einen Integral analytisch zu berechnen.
Mit Methoden der Quadratur können wir Integrale nummerisch berechnen.

\innumpy kann \texttt{scipy.integrate.quad} verwendet werden.
Falls man jedoch eine manuelle Implementation erstellen will, so nutzt man oft die Trapez- oder Simpson-Regel, da sie sowohl einfach zu implementieren, wie auch effizient sind.
In gewissen Anwendungen sind Gauss-Quadratur-Formeln nützlich, welche man durch Spektralmethoden ersetzen kann, welche die FFT verwenden und effizienter sind.

\begin{definition}[]{Quadratur}
    Ein Integral kann durch eine gewichtete Summe von Funktionswerten der Funktion $f$ an verschiedenen Stellen $c_i^n$ approximiert werden:
    \begin{align*}
        \int_{a}^{b} f(x) \dx \approx Q_n(f; a, b) := \sum_{i = 1}^{n} \omega_i^n f(c_i^n)
    \end{align*}
    wobei $\omega_i^n$ die \textit{Gewichte} und $c_i^n \in [a, b]$ die \textit{Knoten} der Quadraturformel sind.
\end{definition}

Wir wollen natürlich wieder $c_i^n \in [a, b]$ und $w_i^n$ so wählen, dass der Fehler minimiert wird.
\begin{definition}[]{Fehler}
    Der Fehler der Quadratur $Q_n(f)$ ist
    \begin{align*}
        E(n) = \left| \int_{a}^{b} f(x) \dx - Q_n(f; a, b) \right|
    \end{align*}
    Wir haben \bi{algebraische Konvergenz} wenn $E(n) = \tco{\frac{1}{n^p}}$ mit $p > 0$ und
    \bi{exponentielle Konvergenz} wenn $E(n) = \tco{q^n}$ mit $0 \leq q < 1$
\end{definition}

% Der Polynom, das Ansatz, yep... excellent German there
Die Idee, den Integral einer schweren Funktion zu berechnen, ist diese mit einer einfachen Funktion, die analytisch integrierbar ist, zu approximieren.
Wenn wir diese Funktion geschickt wählen, dann ist es sogar möglich, dass wir nur eine solche Funktion für alle Funktionen $f$ benötigen.

% NOTE: Called the b_j from the script c_i here, as the clash with the limits of the integral looked a bit confusing
Wir ersetzen also $f$ durch $f_n \in \text{span} \{ c_0, c_1, \ldots, c_n \}$, wobei die $c_i$ eine Basis des Raums der Funktionen auf $[a, b]$ bilden:
\rmvspace
\begin{align*}
    \int_{a}^{b} f(x) \dx x \approx \int_{a}^{b} f_n(x) \dx x = \int_{a}^{b} \left( \sum_{k = 0}^{n} \alpha_k b_k(x) \right) \dx x = \sum_{k = 0}^{n} \alpha_k \int_{a}^{b} c_k(x) \dx x
\end{align*}
Falls wir $c_k(x) = x^k$ haben (was oft der Fall ist, je nach Funktion aber könnte eine rationale Funktion oder andere Arten besser geeignet sein), dann erhalten wir:
\rmvspace
\begin{align*}
    \int_{a}^{b} c_k(x) \dx x = \frac{b^{k + 1} - a^{k + 1}}{k + 1}
\end{align*}
% Mit äquidistanten Punkten $x_0, x_1, \ldots, x_n$ können wir die Koeffizienten $\alpha_0, \alpha_1, \ldots, \alpha_n$ mithilfe eines linearen Gleichunssystems berechnen, 
% wessen Matrix eine $(n + 1) \times (n + 1)$-Vandermonde Matrix ist, so dass $p_n(x) = \alpha_0 + \alpha_1 x + \ldots + \alpha_n x^n$,
% also sind unsere Gleichungen $p_n(x_j) = f(x_j)$ für alle $j = 0, 1, \ldots, n$.
% Dieser Ansatz hat jedoch grosse Rundungsfehler.

\begin{definition}[]{Lagrange-Polynome}
    Für die Knoten $x_0, x_1, \ldots, x_n \in \R$ definieren wir die Polynome
    \rmvspace
    \begin{align*}
        l_i(x) = \prod_{\substack{j = 0\\ j \neq i}}^n \frac{x - x_j}{x_i - x_j}
    \end{align*}

    \rmvspace
    als die Lagrange-Polynome zu den Stützstellen $x_0, x_1, \ldots, x_n$
\end{definition}
Für ein Beispiel verweisen wir auf Beispiel \ref{all:2-3-2}

\setLabelNumber{all}{6}
\fancyremark{Eigenschaften der Lagrange-Polynome} Zu den Eigenschaften aus \ref{all:2-3-4} fügen wir an (die Eigenschaften aus Bemerkung \ref{all:2-3-4} sind hier erneut aufgeführt)
\begin{multicols}{2}
    \begin{enumerate}
        \item $l_i(x_j) = 0 \smallhspace \forall j \neq i$
        \item $l_i(x_i) = 1 \smallhspace \forall i$
        \item $\deg(l_i) = n \smallhspace \forall i$
        \item $\sum_{k = 0}^{n} l_k(x) = 1 \smallhspace \forall x \in \R$
        \item $\sum_{k = 0}^{n} l_k^{(m)}(x) = 0 \text{ für } m > 0$
        \item $l_0, l_1, \ldots, l_n$ bilden Basis von $\mathcal{P}_{n + 1}$
    \end{enumerate}
\end{multicols}

\vspace{-1.5pc}
wobei $\mathcal{P}_{n + 1}$ der Raum der Polynome von Grad maximal $n$ ist.

\fancyremark{Quadraturgewichte aus den Lagrange-Polynomen}
Das Interpolationspolynom ist gegeben durch:
\drmvspace
\begin{align*}
    p(x) = \sum_{j = 0}^{n} f(x_j) l_j(x)
\end{align*}

\drmvspace
Durch die Eigenschaften der Lagrange-Polynome haben wir $p(x_j) = f(x_j)$ und die Konstruktion von $p(x)$ ist eindeutig in $\mathcal{P}_{n + 1}$.
Wir erhalten nun eine Quadraturformel, wenn wir $p$ als Approximation von $f$ verwenden:
\rmvspace
\begin{align*}
    w_j = \int_{a}^{b} l_j(x), \smallhspace j = 0, 1, \ldots, n
\end{align*}

\drmvspace
Durch die Konstruktion der Formel ist sie exakt für alle Polynome aus $\mathcal{P}_{n + 1}$ und der Fehler ist:
\rmvspace
\begin{align*}
    \left| \int_{a}^{b} f(x) \dx x - \int_{a}^{b} p_n(x) \dx x \right| \leq \frac{1}{n!}(b - a)^{n + 1} \max|f^{(n)}(z)|
\end{align*}

\drmvspace
Wir wollen also ein kleines Intervall (oft $b - a < 1$ da wir so das Integral besser approximieren können) und wir setzen voraus, dass $f$ glatt ist.

Da wir aber oft ein grösseres Intervall betrachten möchten, ist ein möglicher Ansatz, das grosse Intervall in kleinere Intervalle zu zerlegen.
Wir nehmen ein äquidistantes Gitter, mit $x_k = x_0 + k \cdot h$ für $h = \frac{b - a}{N}$ und $k = 0, \ldots, N$:
\rmvspace
\begin{align*}
    \int_{a}^{b} f(x) \dx x = \sum_{k = 0}^{N - 1} \int_{x_k}^{x_{k + 1}} f(x) \dx x
\end{align*}

\drmvspace
Die obige Formel wird auch die \textit{summierte} Quadraturformel genannt. Der Fehler ist dann also:
\rmvspace
\begin{align*}
    \left| \int_{a}^{b} f(x) \dx x - \sum_{k = 0}^{N - 1} Q(f, x_k, x_{k + 1}) \right| \leq \ldots \leq C \frac{h^n}{n!}(b - a) &  &
    \text{ mit } C = \max_{z \in [a, b]} |f^{(n)}(z)| = ||f^n||_{\max}
\end{align*}

\rmvspace
Der obige Ansatz ist gewissermassen ``divide and conquer'' (zu Deutsch: ``Teile und Herrsche'', wir werden aber DnC verwenden)
und wir der lokale Fehler liegt in $\tco{h^{n + 1}}$ und mit $N = (b - a) \div h$ Intervallen der Grösse $h$ haben wir einen globalen Fehler in $\tco{h^n}$.
Folglich ist also der Fehler kleiner, je kleiner $h$ ist.

Wir benutzen erneut einen Variablenwechsel, um von einem Referenzintervall $[-1, 1]$ auf eines unserer Teilintervalle $[x_k, x_{k + 1}]$ zu wechseln.
Dies heisst also allgemein für Intervall $[a, b]$ nach $[-1, 1]$:
\begin{align*}
    \int_{a}^{b} f(t) \dx t = \frac{1}{2} (b - a) \int_{-1}^{1} \hat{f}(\tau) \dx \tau & \text{ mit } \hat{f}(\tau) := f\left( \frac{1}{2}(1 - \tau) a + \frac{1}{2}(\tau + 1) b \right)
\end{align*}
Für dieses Referenzintervall können wir die Gewichte $\hat{w}_j$ und die Knoten $\hat{c}_j$ bestimmen.
% OMG, wtf, why can't he decide on using w, \omega or \hat{w} for the weights in the reference interval? That is so dumb.
% FIXME: Choose uniform naming (advocating for \hat{w}_j on the reference intervall and w_j on the [a, b] interval)
\rmvspace
\begin{align*}
    \int_{a}^{b} f(t) \dx t \approx \frac{1}{2}(b - a) \sum_{j = 1}^{n} \hat{w}_j \hat{f}(\hat{c}_j) = \sum_{j = 1}^{n} w_j f(c_j)
     &  & \text{ mit } \elementstack{c_j = \frac{1}{2} (1 - \hat{c}_j) a + \frac{1}{2}(1 + \hat{c}_j) b}{w_j = \frac{1}{2}(b - a)\hat{w}_j}
\end{align*}

\rmvspace\rmvspace
\inlinedef Die Ordnung einer Quadraturformel ist $n$ wenn sie Polynome vom Grad $(n - 1)$ exakt integriert.

Dies folgt natürlich direkt davon, dass wir ein Polynom $n$-ten Grades mit $n + 1$ Koeffizienten darstellen können.

\fancydef{Symmetrie} Eine Quadraturformel auf $[-1, 1]$ heisst symmetrisch, falls $\omega_i = \omega_{n + 1 - i}$ und \\
$c_i = -c_{n + 1 -i}$ gilt für die Gewichte $\omega_i$ und Knoten $c_i$

\inlineremark Die Mittelpunkts-, Trapez- und Simpson-Regeln aus Abschnitt \ref{sec:equidistant-nodes} sind symmetrisch

\inlinetheorem Die Ordnung einer symmetrischen Quadraturformel ist gerade

\inlineproof Kann mittels Induktion bewiesen werden, siehe dazu Seite 123 im Skript