From 3cb2fa354f814fa98474610dac744281285dafc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Patrik=20M=C3=BCller?= Date: Fri, 15 Jul 2022 11:40:55 +0200 Subject: First version of section 'Gauss Quadratur', fix to gamma_approx.py when z=0 --- buch/papers/laguerre/quadratur.tex | 148 ++++++++++++++++++++---- buch/papers/laguerre/references.bib | 11 ++ buch/papers/laguerre/scripts/gamma_approx.py | 18 ++- buch/papers/laguerre/scripts/rel_error_range.py | 2 +- 4 files changed, 151 insertions(+), 28 deletions(-) (limited to 'buch/papers/laguerre') diff --git a/buch/papers/laguerre/quadratur.tex b/buch/papers/laguerre/quadratur.tex index 851fe8a..7cbae48 100644 --- a/buch/papers/laguerre/quadratur.tex +++ b/buch/papers/laguerre/quadratur.tex @@ -5,25 +5,57 @@ % \section{Gauss-Quadratur \label{laguerre:section:quadratur}} - {\large \color{red} TODO: Einleitung und kurze Beschreibung Gauss-Quadratur} - -Siehe Abschnitt~\ref{buch:orthogonalitaet:section:gauss-quadratur} +Die Gauss-Quadratur ist ein numerisches Integrationsverfahren, +welches die Eigenschaften von orthogonalen Polynomen ausnützt. +Herleitungen und Analysen der Gauss-Quadratur können im +Abschnitt~\ref{buch:orthogonalitaet:section:gauss-quadratur} gefunden werden. +Als grundlegende Idee wird die Beobachtung, +dass viele Funktionen sich gut mit Polynomen approximieren lassen, +verwendet. +Stellt man also sicher, +dass ein Verfahren gut für Polynome gut funktioniert, +sollte es auch für andere Funktionen nicht schlecht funktionieren. +Es wird ein Polynom verwendet, +welches an den Punkten $x_0 < x_1 < \ldots < x_n$ +die Funktionwerte~$f(x_i)$ annimmt. +Als Resultat kann das Integral via eine gewichtete Summe der Form \begin{align} \int_a^b f(x) w(x) \, dx \approx \sum_{i=1}^n f(x_i) A_i \label{laguerre:gaussquadratur} \end{align} +berechnet werden. +Die Gauss-Quadratur ist exakt für Polynome mit Grad $2n -1$, +wenn ein Interpolationspolynom von Grad $n$ gewählt wurde. \subsection{Gauss-Laguerre-Quadratur \label{laguerre:subsection:gausslag-quadratur}} -Die Gauss-Quadratur kann auch auf Skalarprodukte mit Gewichtsfunktionen -ausgeweitet werden. -In unserem Falle möchten wir die Gauss Quadratur auf die Laguerre-Polynome -$L_n$ ausweiten. -Diese sind orthogonal im Intervall $(0, \infty)$ bezüglich -der Gewichtsfunktion $e^{-x}$. -Gleichung~\eqref{laguerre:laguerrequadratur} lässt sich wie folgt umformulieren: +Wir möchten nun die Gauss-Quadratur auf die Berechnung +von uneigentlichen Integralen erweitern, +spezifisch auf das Interval $(0, \infty)$. +Mit dem vorher beschriebenen Verfahren ist dies nicht direkt möglich. +Mit einer Transformation die das unendliche Intervall $(a, \infty)$ mit +\begin{align*} +x += +a + \frac{1 - t}{t} +\end{align*} +auf das Intervall $[0, 1]$ transformiert. +Für unser Fall gilt $a = 0$. +Das Integral eines Polynomes in diesem Intervall ist immer divergent, +darum müssen wir sie mit einer Funktion multiplizieren, +die schneller als jedes Polynom gegen $0$ geht, +damit das Integral immer noch konvergiert. +Die Laguerre-Polynome $L_n$ bieten hier Abhilfe, +da ihre Gewichtsfunktion $e^{-x}$ schneller +gegen $0$ konvergiert als jedes Polynom. +% In unserem Falle möchten wir die Gauss Quadratur auf die Laguerre-Polynome +% $L_n$ ausweiten. +% Diese sind orthogonal im Intervall $(0, \infty)$ bezüglich +% der Gewichtsfunktion $e^{-x}$. +Gleichung~\eqref{laguerre:gaussquadratur} lässt sich wie folgt +umformulieren: \begin{align} \int_{0}^{\infty} f(x) e^{-x} dx \approx @@ -43,20 +75,93 @@ l_i(x_j) \delta_{ij} = \begin{cases} -1 & i=j \\ -0 & \text{sonst.} +1 & i=j \\ +0 & \text{.} \end{cases} +% . \end{align*} -Laut \cite{abramowitz+stegun} sind die Gewichte -\begin{align} +die Lagrangschen Interpolationspolynome. +Laut \cite{hildebrand2013introduction} können die Gewicht mit +\begin{align*} A_i + & = +-\frac{C_{n+1} \gamma_n}{C_n \phi'_n(x_i) \phi_{n+1} (x_i)} +\end{align*} +berechnet werden. +$C_i$ entspricht dabei dem Koeffizienten von $x^i$ +des orthogonalen Polynoms $\phi_n(x)$, $\forall i =0,\ldots,n$ und +\begin{align*} +\gamma_n += +\int_0^\infty w(x) \phi_n^2(x)\,dx +\end{align*} +dem Normalisierungsfaktor. +Wir setzen nun $\phi_n(x) = L_n(x)$ und +nutzen den Vorzeichenwechsel der Laguerrekoeffizienten aus, +damit erhalten wir +\begin{align*} +A_i + & = +-\frac{C_{n+1} \gamma_n}{C_n L'_n(x_i) L_{n+1} (x_i)} +\\ + & = \frac{C_n}{C_{n-1}} \frac{\gamma_{n-1}}{L_{n-1}(x_i) L'_n(x_i)} +. +\end{align*} +Für Laguerre-Polynome gilt +\begin{align*} +\frac{C_n}{C_{n-1}} += +-\frac{1}{n} +\quad \text{und} \quad +\gamma_n = -\frac{x_i}{(n + 1)^2 \left[ L_{n + 1}(x_i)\right]^2} +1 +. +\end{align*} +Daraus folgt +\begin{align} +A_i +&= +- \frac{1}{n L_{n-1}(x_i) L'_n(x_i)} +. +\label{laguerre:gewichte_lag_temp} +\end{align} +Nun kann die Rekursionseigenschaft der Laguerre-Polynome +\begin{align*} +x L'_n(x) +&= +n L_n(x) - n L_{n-1}(x) +\\ +&= (x - n - 1) L_n(x) + (n + 1) L_{n+1}(x) +\end{align*} +umgeformt werden und da $x_i$ die Nullstellen von $L_n(x)$ sind, +folgt +\begin{align*} +x_i L'_n(x_i) +&= +- n L_{n-1}(x_i) +\\ +&= + (n + 1) L_{n+1}(x_i) +. +\end{align*} +Setzen wir das nun in \eqref{laguerre:gewichte_lag_temp} ein ergibt sicht +\begin{align} +\nonumber +A_i +&= +\frac{1}{x_i \left[ L'_n(x_i) \right]^2} +\\ +&= +\frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2} . \label{laguerre:quadratur_gewichte} \end{align} \subsubsection{Fehlerterm} +Die Gauss-Laguerre-Quadratur mit $n$ Stützstellen berechnet Integrale +von Polynomen bis zum Grad $2n - 1$ exakt. +Für beliebige Funktionen kann eine Fehlerabschätzung angegeben werden. Der Fehlerterm $R_n$ folgt direkt aus der Approximation \begin{align*} \int_0^{\infty} f(x) e^{-x} \, dx @@ -66,16 +171,15 @@ Der Fehlerterm $R_n$ folgt direkt aus der Approximation und \cite{abramowitz+stegun} gibt ihn als \begin{align} R_n -= + & = +\frac{f^{(2n)}(\xi)}{(2n)!} \int_0^\infty l(x)^2 e^{-x}\,dx +\\ + & = \frac{(n!)^2}{(2n)!} f^{(2n)}(\xi) ,\quad 0 < \xi < \infty \label{laguerre:lag_error} \end{align} an. - -{ -\large \color{red} -TODO: -Noch mehr Text / bessere Beschreibungen in allen Abschnitten -} +Der Fehler ist also abhängig von der $2n$-ten Ableitung +der zu integrierenden Funktion. diff --git a/buch/papers/laguerre/references.bib b/buch/papers/laguerre/references.bib index e12e218..2371922 100644 --- a/buch/papers/laguerre/references.bib +++ b/buch/papers/laguerre/references.bib @@ -4,6 +4,17 @@ % (c) 2020 Autor, Hochschule Rapperswil % +@book{hildebrand2013introduction, + title={Introduction to Numerical Analysis: Second Edition}, + author={Hildebrand, F.B.}, + isbn={9780486318554}, + series={Dover Books on Mathematics}, + url={https://books.google.ch/books?id=ic2jAQAAQBAJ}, + year={2013}, + publisher={Dover Publications}, + pages = {389} +} + @book{abramowitz+stegun, added-at = {2008-06-25T06:25:58.000+0200}, address = {New York}, diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py index 208f770..9f9dae7 100644 --- a/buch/papers/laguerre/scripts/gamma_approx.py +++ b/buch/papers/laguerre/scripts/gamma_approx.py @@ -1,7 +1,5 @@ from pathlib import Path -import matplotlib as mpl -import matplotlib.pyplot as plt import numpy as np import scipy.special @@ -58,6 +56,8 @@ def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11): def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty x, w = _prep_zeros_and_weights(x, w, n) n = len(x) @@ -73,6 +73,8 @@ def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8): def laguerre_gamma_simple(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty x, w = _prep_zeros_and_weights(x, w, n) z += 0j res = np.sum(x ** (z - 1) * w) @@ -81,6 +83,8 @@ def laguerre_gamma_simple(z, x=None, w=None, n=8): def laguerre_gamma_mirror(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty x, w = _prep_zeros_and_weights(x, w, n) z += 0j if z.real < 1e-3: @@ -90,8 +94,8 @@ def laguerre_gamma_mirror(z, x=None, w=None, n=8): return laguerre_gamma_simple(z, x, w) -def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs): - x, w = _prep_zeros_and_weights(x, w, n) +def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs): + x, w = _prep_zeros_and_weights(x, w, n) if func == "simple": f = laguerre_gamma_simple elif func == "mirror": @@ -104,4 +108,8 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs): def calc_rel_error(x, y): - return (y - x) / x + mask = np.abs(x) != np.infty + rel_error = np.zeros_like(y) + rel_error[mask] = (y[mask] - x[mask]) / x[mask] + rel_error[~mask] = 0.0 + return rel_error diff --git a/buch/papers/laguerre/scripts/rel_error_range.py b/buch/papers/laguerre/scripts/rel_error_range.py index 7d017a7..7c74d76 100644 --- a/buch/papers/laguerre/scripts/rel_error_range.py +++ b/buch/papers/laguerre/scripts/rel_error_range.py @@ -5,7 +5,7 @@ if __name__ == "__main__": import gamma_approx as ga - N = 1000 + N = 1001 xmin = -5 xmax = 5 ns = np.arange(2, 12, 2) -- cgit v1.2.1