aboutsummaryrefslogtreecommitdiffstats
path: root/buch/papers/laguerre
diff options
context:
space:
mode:
Diffstat (limited to 'buch/papers/laguerre')
-rw-r--r--buch/papers/laguerre/Makefile6
-rw-r--r--buch/papers/laguerre/Makefile.inc5
-rw-r--r--buch/papers/laguerre/definition.tex177
-rw-r--r--buch/papers/laguerre/eigenschaften.tex115
-rw-r--r--buch/papers/laguerre/gamma.tex76
-rw-r--r--buch/papers/laguerre/images/laguerre_polynomes.pdfbin0 -> 16239 bytes
-rw-r--r--buch/papers/laguerre/main.tex7
-rw-r--r--buch/papers/laguerre/packages.tex1
-rw-r--r--buch/papers/laguerre/quadratur.tex78
-rw-r--r--buch/papers/laguerre/references.bib45
-rw-r--r--buch/papers/laguerre/scripts/gamma_approx.ipynb395
-rw-r--r--buch/papers/laguerre/scripts/laguerre_plot.py100
-rw-r--r--buch/papers/laguerre/transformation.tex31
-rw-r--r--buch/papers/laguerre/wasserstoff.tex29
14 files changed, 917 insertions, 148 deletions
diff --git a/buch/papers/laguerre/Makefile b/buch/papers/laguerre/Makefile
index 606d7e1..0f0985a 100644
--- a/buch/papers/laguerre/Makefile
+++ b/buch/papers/laguerre/Makefile
@@ -4,6 +4,8 @@
# (c) 2020 Prof Dr Andreas Mueller
#
-images:
- @echo "no images to be created in laguerre"
+images: images/laguerre_polynomes.pdf
+
+images/laguerre_polynomes.pdf: scripts/laguerre_plot.py
+ python3 scripts/laguerre_plot.py
diff --git a/buch/papers/laguerre/Makefile.inc b/buch/papers/laguerre/Makefile.inc
index 1eb5034..12b0935 100644
--- a/buch/papers/laguerre/Makefile.inc
+++ b/buch/papers/laguerre/Makefile.inc
@@ -9,8 +9,7 @@ dependencies-laguerre = \
papers/laguerre/references.bib \
papers/laguerre/definition.tex \
papers/laguerre/eigenschaften.tex \
- papers/laguerre/quadratur.tex \
- papers/laguerre/transformation.tex \
- papers/laguerre/wasserstoff.tex
+ papers/laguerre/quadratur.tex \
+ papers/laguerre/gamma.tex
diff --git a/buch/papers/laguerre/definition.tex b/buch/papers/laguerre/definition.tex
index 5f6d8bd..d111f6f 100644
--- a/buch/papers/laguerre/definition.tex
+++ b/buch/papers/laguerre/definition.tex
@@ -4,45 +4,154 @@
% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
%
\section{Definition
-\label{laguerre:section:definition}}
+ \label{laguerre:section:definition}}
\rhead{Definition}
-
+Die verallgemeinerte Laguerre-Differentialgleichung ist gegeben durch
\begin{align}
- x y''(x) + (1 - x) y'(x) + n y(x)
- =
- 0
- \label{laguerre:dgl}
+x y''(x) + (\nu + 1 - x) y'(x) + n y(x)
+=
+0
+, \quad
+n \in \mathbb{N}_0
+, \quad
+x \in \mathbb{R}
+.
+\label{laguerre:dgl}
\end{align}
-
+Die klassische Laguerre-Diffentialgleichung erhält man, wenn $\nu = 0$.
+Hier wird die verallgemeinerte Laguerre-Differentialgleichung verwendet,
+weil die Lösung mit der selben Methode berechnet werden kann,
+aber man zusätzlich die Lösung für den allgmeinen Fall erhält.
+Zur Lösung der Gleichung \eqref{laguerre:dgl} verwenden wir einen
+Potenzreihenansatz.
+Da wir bereits wissen, dass die Lösung orthogonale Polynome sind,
+erscheint dieser Ansatz sinnvoll.
+Setzt man nun den Ansatz
+\begin{align*}
+y(x)
+ & =
+\sum_{k=0}^\infty a_k x^k
+\\
+y'(x)
+ & =
+\sum_{k=1}^\infty k a_k x^{k-1}
+=
+\sum_{k=0}^\infty (k+1) a_{k+1} x^k
+\\
+y''(x)
+ & =
+\sum_{k=2}^\infty k (k-1) a_k x^{k-2}
+=
+\sum_{k=1}^\infty (k+1) k a_{k+1} x^{k-1}
+\end{align*}
+in die Differentialgleichung ein, erhält man:
+\begin{align*}
+\sum_{k=1}^\infty (k+1) k a_{k+1} x^k
++
+(\nu + 1)\sum_{k=0}^\infty (k+1) a_{k+1} x^k
+-
+\sum_{k=0}^\infty k a_k x^k
++
+n \sum_{k=0}^\infty a_k x^k
+ & =
+0 \\
+\sum_{k=1}^\infty
+\left[ (k+1) k a_{k+1} + (\nu + 1)(k+1) a_{k+1} - k a_k + n a_k \right] x^k
+ & =
+0.
+\end{align*}
+Daraus lässt sich die Rekursionsbeziehung
+\begin{align*}
+a_{k+1}
+ & =
+\frac{k-n}{(k+1) (k + \nu + 1)} a_k
+\end{align*}
+ableiten.
+Für ein konstantes $n$ erhalten wir als Potenzreihenlösung ein Polynom vom Grad
+$n$,
+denn für $k=n$ wird $a_{n+1} = 0$ und damit auch $a_{n+2}=a_{n+3}=\ldots=0$.
+Aus der Rekursionsbeziehung ist zudem ersichtlich,
+dass $a_0 \neq 0$ beliebig gewählt werden kann.
+Wählen wir nun $a_0 = 1$, dann folgt für die Koeffizienten $a_1, a_2, a_3$
+\begin{align*}
+a_1
+=
+-\frac{n}{1 \cdot (\nu + 1)}
+, & &
+a_2
+=
+\frac{(n-1)n}{1 \cdot 2 \cdot (\nu + 1)(\nu + 2)}
+, & &
+a_3
+=
+-\frac{(n-2)(n-1)n}{1 \cdot 2 \cdot 3 \cdot (\nu + 1)(\nu + 2)(\nu + 3)}
+\end{align*}
+und allgemein
+\begin{align*}
+k
+ & \leq
+n:
+ &
+a_k
+ & =
+(-1)^k \frac{n!}{(n-k)!} \frac{1}{k!(\nu + 1)_k}
+=
+\frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k}
+\\
+k & >n:
+ &
+a_k
+ & =
+0.
+\end{align*}
+Somit erhalten wir für $\nu = 0$ die Laguerre-Polynome
\begin{align}
- L_n(x)
- =
- \sum_{k=0}^{n}
- \frac{(-1)^k}{k!}
- \begin{pmatrix}
- n \\
- k
- \end{pmatrix}
- x^k
- \label{laguerre:polynom}
+L_n(x)
+=
+\sum_{k=0}^{n} \frac{(-1)^k}{k!} \binom{n}{k} x^k
+\label{laguerre:polynom}
\end{align}
-
+und mit $\nu \in \mathbb{R}$ die verallgemeinerten Laguerre-Polynome
\begin{align}
- x y''(x) + (\alpha + 1 - x) y'(x) + n y(x)
- =
- 0
- \label{laguerre:generell_dgl}
+L_n^\nu(x)
+=
+\sum_{k=0}^{n} \frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k} x^k.
+\label{laguerre:allg_polynom}
\end{align}
-\begin{align}
- L_n^\alpha (x)
- =
- \sum_{k=0}^{n}
- \frac{(-1)^k}{k!}
- \begin{pmatrix}
- n + \alpha \\
- n - k
- \end{pmatrix}
- x^k
- \label{laguerre:polynom}
-\end{align}
+\subsection{Analytische Fortsetzung}
+Durch die analytische Fortsetzung erhalten wir zudem noch die zweite Lösung der
+Differentialgleichung mit der Form
+\begin{align*}
+\Xi_n(x)
+=
+L_n(x) \ln(x) + \sum_{k=1}^\infty d_k x^k
+\end{align*}
+Nach einigen mühsamen Rechnungen,
+die den Rahmen dieses Kapitel sprengen würden,
+erhalten wir
+\begin{align*}
+\Xi_n
+=
+L_n(x) \ln(x)
++
+\sum_{k=1}^n \frac{(-1)^k}{k!} \binom{n}{k}
+(\alpha_{n-k} - \alpha_n - 2 \alpha_k)x^k
++
+(-1)^n \sum_{k=1}^\infty \frac{(k-1)!n!}{((n+k)!)^2} x^{n+k},
+\end{align*}
+wobei $\alpha_0 = 0$ und $\alpha_k =\sum_{i=1}^k i^{-1}$,
+$\forall k \in \mathbb{N}$.
+Die Laguerre-Polynome von Grad $0$ bis $7$ sind in
+Abbildung~\ref{laguerre:fig:polyeval} dargestellt.
+\begin{figure}
+\centering
+\includegraphics[width=0.7\textwidth]{%
+ papers/laguerre/images/laguerre_polynomes.pdf%
+}
+\caption{Laguerre-Polynome vom Grad $0$ bis $7$}
+\label{laguerre:fig:polyeval}
+\end{figure}
+
+% https://www.math.kit.edu/iana1/lehre/hm3phys2012w/media/laguerre.pdf
+% http://www.physics.okayama-u.ac.jp/jeschke_homepage/E4/kapitel4.pdf
diff --git a/buch/papers/laguerre/eigenschaften.tex b/buch/papers/laguerre/eigenschaften.tex
index b7597e5..b0cc3a3 100644
--- a/buch/papers/laguerre/eigenschaften.tex
+++ b/buch/papers/laguerre/eigenschaften.tex
@@ -4,5 +4,116 @@
% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
%
\section{Eigenschaften
-\label{laguerre:section:eigenschaften}}
-\rhead{Eigenschaften} \ No newline at end of file
+ \label{laguerre:section:eigenschaften}}
+{
+\large \color{red}
+TODO:
+Evtl. nur Orthogonalität hier behandeln, da nur diese für die Gauss-Quadratur
+benötigt wird.
+}
+
+Die Laguerre-Polynome besitzen einige interessante Eigenschaften
+\rhead{Eigenschaften}
+
+\subsection{Orthogonalität
+ \label{laguerre:subsection:orthogonal}}
+Im Abschnitt~\ref{laguerre:section:definition} haben wir behauptet,
+dass die Laguerre-Polynome orthogonale Polynome sind.
+Zu dieser Behauptung möchten wir nun einen Beweis liefern.
+Wenn wir die Laguerre\--Differentialgleichung in ein
+Sturm\--Liouville\--Problem umwandeln können, haben wir bewiesen, dass es sich
+bei
+den Laguerre\--Polynomen um orthogonale Polynome handelt (siehe
+Abschnitt~\ref{buch:integrale:subsection:sturm-liouville-problem}).
+Der Sturm-Liouville-Operator hat die Form
+\begin{align}
+S
+=
+\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right).
+\label{laguerre:slop}
+\end{align}
+Aus der Beziehung
+\begin{align}
+S
+ & =
+\Lambda
+\nonumber
+\\
+\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right)
+ & =
+x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx}
+\label{laguerre:sl-lag}
+\end{align}
+lässt sich sofort erkennen, dass $q(x) = 0$.
+Ausserdem ist ersichtlich, dass $p(x)$ die Differentialgleichung
+\begin{align*}
+x \frac{dp}{dx}
+=
+-(\nu + 1 - x) p,
+\end{align*}
+erfüllen muss.
+Durch Separation erhalten wir dann
+\begin{align*}
+\int \frac{dp}{p}
+ & =
+-\int \frac{\nu + 1 - x}{x}dx
+\\
+\log p
+ & =
+-\log \nu + 1 - x + C
+\\
+p(x)
+ & =
+-C x^{\nu + 1} e^{-x}
+\end{align*}
+Eingefügt in Gleichung~\eqref{laguerre:sl-lag} erhalten wir
+\begin{align*}
+\frac{C}{w(x)}
+\left(
+x^{\nu+1} e^{-x} \frac{d^2}{dx^2} +
+(\nu + 1 - x) x^{\nu} e^{-x} \frac{d}{dx}
+\right)
+=
+x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx}.
+\end{align*}
+Mittels Koeffizientenvergleich kann nun abgelesen werden, dass $w(x) = x^\nu
+e^{-x}$ und $C=1$ mit $\nu > -1$.
+Die Gewichtsfunktion $w(x)$ wächst für $x\rightarrow-\infty$ sehr schnell an,
+deshalb ist die Laguerre-Gewichtsfunktion nur geeignet für den
+Definitionsbereich $(0, \infty)$.
+Bleibt nur noch sicherzustellen, dass die Randbedingungen,
+\begin{align}
+k_0 y(0) + h_0 p(0)y'(0)
+ & =
+0
+\label{laguerre:sllag_randa}
+\\
+k_\infty y(\infty) + h_\infty p(\infty) y'(\infty)
+ & =
+0
+\label{laguerre:sllag_randb}
+\end{align}
+mit $|k_i|^2 + |h_i|^2 \neq 0,\,\forall i \in \{0, \infty\}$, erfüllt sind.
+Am linken Rand (Gleichung~\eqref{laguerre:sllag_randa}) kann $y(0) = 1$, $k_0 =
+0$ und $h_0 = 1$ verwendet werden,
+was auch die Laguerre-Polynome ergeben haben.
+Für den rechten Rand ist die Bedingung (Gleichung~\eqref{laguerre:sllag_randb})
+\begin{align*}
+\lim_{x \rightarrow \infty} p(x) y'(x)
+ & =
+\lim_{x \rightarrow \infty} -x^{\nu + 1} e^{-x} y'(x)
+=
+0
+\end{align*}
+für beliebige Polynomlösungen erfüllt für $k_\infty=0$ und $h_\infty=1$.
+Damit können wir schlussfolgern, dass die Laguerre-Polynome orthogonal
+bezüglich des Skalarproduktes auf dem Intervall $(0, \infty)$ mit der Laguerre\--Gewichtsfunktion
+$w(x)=x^\nu e^{-x}$ sind.
+
+
+\subsection{Rodrigues-Formel}
+
+\subsection{Drei-Terme Rekursion}
+
+\subsection{Beziehung mit der Hypergeometrischen Funktion}
+
diff --git a/buch/papers/laguerre/gamma.tex b/buch/papers/laguerre/gamma.tex
new file mode 100644
index 0000000..e3838b0
--- /dev/null
+++ b/buch/papers/laguerre/gamma.tex
@@ -0,0 +1,76 @@
+%
+% gamma.tex
+%
+% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
+%
+\section{Anwendung: Berechnung der Gamma-Funktion
+ \label{laguerre:section:quad-gamma}}
+Die Gauss-Laguerre-Quadratur kann nun verwendet werden,
+um exponentiell abfallende Funktionen im Definitionsbereich $(0, \infty)$ zu
+berechnen.
+Dabei bietet sich z.B. die Gamma-Funkion bestens an, wie wir in den folgenden
+Abschnitten sehen werden.
+
+\subsection{Gamma-Funktion}
+Die Gamma-Funktion ist eine Erweiterung der Fakultät auf die reale und komplexe
+Zahlenmenge.
+Die Definition~\ref{buch:rekursion:def:gamma} beschreibt die Gamma-Funktion als
+Integral der Form
+\begin{align}
+\Gamma(z)
+ & =
+\int_0^\infty t^{z-1} e^{-t} dt
+,
+\quad
+\text{wobei Realteil von $z$ grösser als $0$}
+,
+\label{laguerre:gamma}
+\end{align}
+welches alle Eigenschaften erfüllt, um mit der Gauss-Laguerre-Quadratur
+berechnet zu werden.
+
+\subsubsection{Funktionalgleichung}
+Die Funktionalgleichung besagt
+\begin{align}
+z \Gamma(z) = \Gamma(z+1).
+\label{laguerre:gamma_funktional}
+\end{align}
+Mittels dieser Gleichung kann der Wert an einer bestimmten,
+geeigneten Stelle evaluiert werden und dann zurückverschoben werden,
+um das gewünschte Resultat zu erhalten.
+
+\subsection{Berechnung mittels Gauss-Laguerre-Quadratur}
+
+Fehlerterm:
+\begin{align*}
+R_n
+=
+(z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z-2n-1}
+\end{align*}
+
+\subsubsection{Finden der optimalen Berechnungsstelle}
+Nun stellt sich die Frage,
+ob die Approximation mittels Gauss-Laguerre-Quadratur verbessert werden kann,
+wenn man das Problem an einer geeigneten Stelle evaluiert und
+dann zurückverschiebt mit der Funktionalgleichung.
+Dazu wollen wir den Fehlerterm in
+Gleichung~\eqref{laguerre:lagurre:lag_error} anpassen und dann minimieren.
+Zunächst wollen wir dies nur für $z\in \mathbb{R}$ und $0<z<1$ definieren.
+Zudem nehmen wir an, dass die optimale Stelle $x^* \in \mathbb{R}$, $z < x^*$
+ist.
+Dann fügen wir einen Verschiebungsterm um $m$ Stellen ein, daraus folgt
+\begin{align*}
+R_n
+=
+\frac{(z - 2n)_{2n}}{(z - m)_m} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1}
+.
+\end{align*}
+
+{
+\large \color{red}
+TODO:
+Geeignete Minimierung für Fehler finden, so dass sie mit den emprisich
+bestimmen optimalen Punkten übereinstimmen.
+}
+
+\subsection{Resultate}
diff --git a/buch/papers/laguerre/images/laguerre_polynomes.pdf b/buch/papers/laguerre/images/laguerre_polynomes.pdf
new file mode 100644
index 0000000..3976bc7
--- /dev/null
+++ b/buch/papers/laguerre/images/laguerre_polynomes.pdf
Binary files differ
diff --git a/buch/papers/laguerre/main.tex b/buch/papers/laguerre/main.tex
index 1fe0f8b..00e3b43 100644
--- a/buch/papers/laguerre/main.tex
+++ b/buch/papers/laguerre/main.tex
@@ -8,13 +8,14 @@
\begin{refsection}
\chapterauthor{Patrik Müller}
-Hier kommt eine Einleitung.
+{\large \color{red} TODO: Einleitung}
\input{papers/laguerre/definition}
\input{papers/laguerre/eigenschaften}
\input{papers/laguerre/quadratur}
-\input{papers/laguerre/transformation}
-\input{papers/laguerre/wasserstoff}
+\input{papers/laguerre/gamma}
+% \input{papers/laguerre/transformation}
+% \input{papers/laguerre/wasserstoff}
\printbibliography[heading=subbibliography]
\end{refsection}
diff --git a/buch/papers/laguerre/packages.tex b/buch/papers/laguerre/packages.tex
index ab55228..4ebc172 100644
--- a/buch/papers/laguerre/packages.tex
+++ b/buch/papers/laguerre/packages.tex
@@ -7,4 +7,3 @@
% if your paper needs special packages, add package commands as in the
% following example
\usepackage{derivative}
-
diff --git a/buch/papers/laguerre/quadratur.tex b/buch/papers/laguerre/quadratur.tex
index 8ab1af5..60fad7f 100644
--- a/buch/papers/laguerre/quadratur.tex
+++ b/buch/papers/laguerre/quadratur.tex
@@ -3,27 +3,77 @@
%
% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
%
-\section{Gauss-Laguerre Quadratur
-\label{laguerre:section:quadratur}}
+\section{Gauss-Quadratur
+ \label{laguerre:section:quadratur}}
+ {\large \color{red} TODO: Einleitung und kurze Beschreibung Gauss-Quadratur}
+\begin{align}
+\int_a^b f(x) w(x)
+\approx
+\sum_{i=1}^N f(x_i) A_i
+\label{laguerre:gaussquadratur}
+\end{align}
+\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 wiefolgt umformulieren:
\begin{align}
- \int_a^b f(x) w(x)
- \approx
- \sum_{i=1}^N f(x_i) A_i
- \label{laguerre:gaussquadratur}
+\int_{0}^{\infty} f(x) e^{-x} dx
+\approx
+\sum_{i=1}^{N} f(x_i) A_i
+\label{laguerre:laguerrequadratur}
\end{align}
+\subsubsection{Stützstellen und Gewichte}
+Nach der Definition der Gauss-Quadratur müssen als Stützstellen die Nullstellen
+des verwendeten Polynoms genommen werden.
+Das heisst für das Laguerre-Polynom $L_n$ müssen dessen Nullstellen $x_i$ und
+als Gewichte $A_i$ werden die Integrale $l_i(x)e^{-x}$ verwendet werden.
+Dabei sind
+\begin{align*}
+l_i(x_j)
+=
+\delta_{ij}
+=
+\begin{cases}
+1 & i=j \\
+0 & \text{sonst.}
+\end{cases}
+\end{align*}
+Laut \cite{abramowitz+stegun} sind die Gewichte also
\begin{align}
- \int_{0}^{\infty} f(x) e^{-x} dx
- \approx
- \sum_{i=1}^{N} f(x_i) A_i
- \label{laguerre:laguerrequadratur}
+A_i
+=
+\frac{x_i}{(n + 1)^2 \left[ L_{n + 1}(x_i)\right]^2}
+.
+\label{laguerre:quadratur_gewichte}
\end{align}
+\subsubsection{Fehlerterm}
+Der Fehlerterm $R_n$ folgt direkt aus der Approximation
+\begin{align*}
+\int_0^{\infty} f(x) e^{-x} dx
+=
+\sum_{i=1}^n f(x_i) A_i + R_n
+\end{align*}
+un \cite{abramowitz+stegun} gibt in als
\begin{align}
- A_i
- =
- \frac{x_i}{(n + 1)^2 \left[ L_{n + 1}(x_i)\right]^2}
- \label{laguerre:quadratur_gewichte}
+R_n
+=
+\frac{(n!)^2}{(2n)!} f^{(2n)}(\xi)
+,\quad
+0 < \xi < \infty
+\label{lagurre:lag_error}
\end{align}
+an.
+{
+\large \color{red}
+TODO:
+Noch mehr Text / bessere Beschreibungen in allen Abschnitten
+}
diff --git a/buch/papers/laguerre/references.bib b/buch/papers/laguerre/references.bib
index caf270f..6956ade 100644
--- a/buch/papers/laguerre/references.bib
+++ b/buch/papers/laguerre/references.bib
@@ -4,32 +4,19 @@
% (c) 2020 Autor, Hochschule Rapperswil
%
-@online{laguerre:bibtex,
- title = {BibTeX},
- url = {https://de.wikipedia.org/wiki/BibTeX},
- date = {2020-02-06},
- year = {2020},
- month = {2},
- day = {6}
-}
-
-@book{laguerre:numerical-analysis,
- title = {Numerical Analysis},
- author = {David Kincaid and Ward Cheney},
- publisher = {American Mathematical Society},
- year = {2002},
- isbn = {978-8-8218-4788-6},
- inseries = {Pure and applied undegraduate texts},
- volume = {2}
-}
-
-@article{laguerre:mendezmueller,
- author = { Tabea Méndez and Andreas Müller },
- title = { Noncommutative harmonic analysis and image registration },
- journal = { Appl. Comput. Harmon. Anal.},
- year = 2019,
- volume = 47,
- pages = {607--627},
- url = {https://doi.org/10.1016/j.acha.2017.11.004}
-}
-
+@book{abramowitz+stegun,
+ added-at = {2008-06-25T06:25:58.000+0200},
+ address = {New York},
+ author = {Abramowitz, Milton and Stegun, Irene A.},
+ biburl = {https://www.bibsonomy.org/bibtex/223ec744709b3a776a1af0a3fd65cd09f/a_olympia},
+ description = {BibTeX - Wikipedia, the free encyclopedia},
+ edition = {ninth Dover printing, tenth GPO printing},
+ interhash = {d4914a420f489f7c5129ed01ec3cf80c},
+ intrahash = {23ec744709b3a776a1af0a3fd65cd09f},
+ keywords = {Handbook},
+ publisher = {Dover},
+ pages = {890},
+ timestamp = {2008-06-25T06:25:58.000+0200},
+ title = {Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables},
+ year = 1972
+} \ No newline at end of file
diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb
new file mode 100644
index 0000000..44f3abd
--- /dev/null
+++ b/buch/papers/laguerre/scripts/gamma_approx.ipynb
@@ -0,0 +1,395 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Gauss-Laguerre Quadratur für die Gamma-Funktion\n",
+ "\n",
+ "$$\n",
+ " \\Gamma(z)\n",
+ " = \n",
+ " \\int_0^\\infty t^{z-1}e^{-t}dt\n",
+ "$$\n",
+ "\n",
+ "$$\n",
+ " \\int_0^\\infty f(x) e^{-x} dx \n",
+ " \\approx \n",
+ " \\sum_{i=1}^{N} f(x_i) w_i\n",
+ " \\qquad\\text{ wobei }\n",
+ " w_i = \\frac{x_i}{(n+1)^2 [L_{n+1}(x_i)]^2}\n",
+ "$$\n",
+ "und $x_i$ sind Nullstellen des Laguerre Polynoms $L_n(x)$\n",
+ "\n",
+ "Der Fehler ist gegeben als\n",
+ "\n",
+ "$$\n",
+ " E \n",
+ " =\n",
+ " \\frac{(n!)^2}{(2n)!} f^{(2n)}(\\xi) \n",
+ " = \n",
+ " \\frac{(-2n + z)_{2n}}{(z-m)_m} \\frac{(n!)^2}{(2n)!} \\xi^{z + m - 2n - 1}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from cmath import exp, pi, sin, sqrt\n",
+ "import scipy.special\n",
+ "\n",
+ "EPSILON = 1e-07\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "lanczos_p = [\n",
+ " 676.5203681218851,\n",
+ " -1259.1392167224028,\n",
+ " 771.32342877765313,\n",
+ " -176.61502916214059,\n",
+ " 12.507343278686905,\n",
+ " -0.13857109526572012,\n",
+ " 9.9843695780195716e-6,\n",
+ " 1.5056327351493116e-7,\n",
+ "]\n",
+ "\n",
+ "\n",
+ "def drop_imag(z):\n",
+ " if abs(z.imag) <= EPSILON:\n",
+ " z = z.real\n",
+ " return z\n",
+ "\n",
+ "\n",
+ "def lanczos_gamma(z):\n",
+ " z = complex(z)\n",
+ " if z.real < 0.5:\n",
+ " y = pi / (sin(pi * z) * lanczos_gamma(1 - z)) # Reflection formula\n",
+ " else:\n",
+ " z -= 1\n",
+ " x = 0.99999999999980993\n",
+ " for (i, pval) in enumerate(lanczos_p):\n",
+ " x += pval / (z + i + 1)\n",
+ " t = z + len(lanczos_p) - 0.5\n",
+ " y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x\n",
+ " return drop_imag(y)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "zeros, weights = np.polynomial.laguerre.laggauss(8)\n",
+ "# zeros = np.array(\n",
+ "# [\n",
+ "# 1.70279632305101000e-1,\n",
+ "# 9.03701776799379912e-1,\n",
+ "# 2.25108662986613069e0,\n",
+ "# 4.26670017028765879e0,\n",
+ "# 7.04590540239346570e0,\n",
+ "# 1.07585160101809952e1,\n",
+ "# 1.57406786412780046e1,\n",
+ "# 2.28631317368892641e1,\n",
+ "# ]\n",
+ "# )\n",
+ "\n",
+ "# weights = np.array(\n",
+ "# [\n",
+ "# 3.69188589341637530e-1,\n",
+ "# 4.18786780814342956e-1,\n",
+ "# 1.75794986637171806e-1,\n",
+ "# 3.33434922612156515e-2,\n",
+ "# 2.79453623522567252e-3,\n",
+ "# 9.07650877335821310e-5,\n",
+ "# 8.48574671627253154e-7,\n",
+ "# 1.04800117487151038e-9,\n",
+ "# ]\n",
+ "# )\n",
+ "\n",
+ "\n",
+ "def pochhammer(z, n):\n",
+ " return np.prod(z + np.arange(n))\n",
+ "\n",
+ "\n",
+ "def find_shift(z, target):\n",
+ " factor = 1.0\n",
+ " steps = int(np.floor(target - np.real(z)))\n",
+ " zs = z + steps\n",
+ " if steps > 0:\n",
+ " factor = 1 / pochhammer(z, steps)\n",
+ " elif steps < 0:\n",
+ " factor = pochhammer(zs, -steps)\n",
+ " return zs, factor\n",
+ "\n",
+ "\n",
+ "def laguerre_gamma(z, x, w, target=11):\n",
+ " # res = 0.0\n",
+ " z = complex(z)\n",
+ " if z.real < 1e-3:\n",
+ " res = pi / (\n",
+ " sin(pi * z) * laguerre_gamma(1 - z, x, w, target)\n",
+ " ) # Reflection formula\n",
+ " else:\n",
+ " z_shifted, correction_factor = find_shift(z, target)\n",
+ " res = np.sum(x ** (z_shifted - 1) * w)\n",
+ " res *= correction_factor\n",
+ " res = drop_imag(res)\n",
+ " return res\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def eval_laguerre(x, target=12):\n",
+ " return np.array([laguerre_gamma(xi, zeros, weights, target) for xi in x])\n",
+ "\n",
+ "\n",
+ "def eval_lanczos(x):\n",
+ " return np.array([lanczos_gamma(xi) for xi in x])\n",
+ "\n",
+ "\n",
+ "def eval_mean_laguerre(x, targets):\n",
+ " return np.mean([eval_laguerre(x, target) for target in targets], 0)\n",
+ "\n",
+ "\n",
+ "def calc_rel_error(x, y):\n",
+ " return (y - x) / x\n",
+ "\n",
+ "\n",
+ "def evaluate(x, target=12):\n",
+ " lanczos_gammas = eval_lanczos(x)\n",
+ " laguerre_gammas = eval_laguerre(x, target)\n",
+ " rel_error = calc_rel_error(lanczos_gammas, laguerre_gammas)\n",
+ " return rel_error\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Test with real values"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Empirische Tests zeigen:\n",
+ "- $n=4 \\Rightarrow m=6$\n",
+ "- $n=5 \\Rightarrow m=7$ oder $m=8$\n",
+ "- $n=6 \\Rightarrow m=9$\n",
+ "- $n=7 \\Rightarrow m=10$\n",
+ "- $n=8 \\Rightarrow m=11$ oder $m=12$\n",
+ "- $n=9 \\Rightarrow m=13$\n",
+ "- $n=10 \\Rightarrow m=14$\n",
+ "- $n=11 \\Rightarrow m=15$ oder $m=16$\n",
+ "- $n=12 \\Rightarrow m=17$\n",
+ "- $n=13 \\Rightarrow m=18 \\Rightarrow $ Beginnt numerisch instabil zu werden \n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "zeros, weights = np.polynomial.laguerre.laggauss(12)\n",
+ "targets = np.arange(16, 21)\n",
+ "mean_targets = ((16, 17),)\n",
+ "x = np.linspace(EPSILON, 1 - EPSILON, 101)\n",
+ "_, axs = plt.subplots(\n",
+ " 2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n",
+ ")\n",
+ "\n",
+ "lanczos = eval_lanczos(x)\n",
+ "for mean_target in mean_targets:\n",
+ " vals = eval_mean_laguerre(x, mean_target)\n",
+ " rel_error_mean = calc_rel_error(lanczos, vals)\n",
+ " axs[0].plot(x, rel_error_mean, label=mean_target)\n",
+ " axs[1].semilogy(x, np.abs(rel_error_mean), label=mean_target)\n",
+ "\n",
+ "mins = []\n",
+ "maxs = []\n",
+ "for target in targets:\n",
+ " rel_error = evaluate(x, target)\n",
+ " mins.append(np.min(np.abs(rel_error[(0.1 <= x) & (x <= 0.9)])))\n",
+ " maxs.append(np.max(np.abs(rel_error)))\n",
+ " axs[0].plot(x, rel_error, label=target)\n",
+ " axs[1].semilogy(x, np.abs(rel_error), label=target)\n",
+ "# axs[0].set_ylim(*(np.array([-1, 1]) * 3.5e-8))\n",
+ "\n",
+ "axs[0].set_xlim(x[0], x[-1])\n",
+ "axs[1].set_ylim(np.min(mins), 1.04*np.max(maxs))\n",
+ "for ax in axs:\n",
+ " ax.legend()\n",
+ " ax.grid(which=\"both\")\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "targets = (16, 17)\n",
+ "xmax = 15\n",
+ "x = np.linspace(-xmax + EPSILON, xmax - EPSILON, 1000)\n",
+ "\n",
+ "mean_lag = eval_mean_laguerre(x, targets)\n",
+ "lanczos = eval_lanczos(x)\n",
+ "rel_error = calc_rel_error(lanczos, mean_lag)\n",
+ "rel_error_simple = evaluate(x, targets[-1])\n",
+ "# rel_error = evaluate(x, target)\n",
+ "\n",
+ "_, axs = plt.subplots(\n",
+ " 2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n",
+ ")\n",
+ "axs[0].plot(x, rel_error, label=targets)\n",
+ "axs[1].semilogy(x, np.abs(rel_error), label=targets)\n",
+ "axs[0].plot(x, rel_error_simple, label=targets[-1])\n",
+ "axs[1].semilogy(x, np.abs(rel_error_simple), label=targets[-1])\n",
+ "axs[0].set_xlim(x[0], x[-1])\n",
+ "# axs[0].set_ylim(*(np.array([-1, 1]) * 4.2e-8))\n",
+ "# axs[1].set_ylim(1e-10, 5e-8)\n",
+ "for ax in axs:\n",
+ " ax.legend()\n",
+ "\n",
+ "x2 = np.linspace(-5 + EPSILON, 5, 4001)\n",
+ "_, ax = plt.subplots(constrained_layout=True, figsize=(8, 6))\n",
+ "ax.plot(x2, eval_mean_laguerre(x2, targets))\n",
+ "ax.set_xlim(x2[0], x2[-1])\n",
+ "ax.set_ylim(-7.5, 25)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Test with complex values"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "targets = (16, 17)\n",
+ "vals = np.linspace(-5 + EPSILON, 5, 100)\n",
+ "x, y = np.meshgrid(vals, vals)\n",
+ "mesh = x + 1j * y\n",
+ "input = mesh.flatten()\n",
+ "\n",
+ "mean_lag = eval_mean_laguerre(input, targets).reshape(mesh.shape)\n",
+ "lanczos = eval_lanczos(input).reshape(mesh.shape)\n",
+ "rel_error = np.abs(calc_rel_error(lanczos, mean_lag))\n",
+ "\n",
+ "lag = eval_laguerre(input, targets[-1]).reshape(mesh.shape)\n",
+ "rel_error_simple = np.abs(calc_rel_error(lanczos, lag))\n",
+ "# rel_error = evaluate(x, target)\n",
+ "\n",
+ "fig, axs = plt.subplots(\n",
+ " 2,\n",
+ " 2,\n",
+ " sharex=True,\n",
+ " sharey=True,\n",
+ " clear=True,\n",
+ " constrained_layout=True,\n",
+ " figsize=(12, 10),\n",
+ ")\n",
+ "_c = axs[0, 1].pcolormesh(x, y, np.log10(np.abs(lanczos - mean_lag)), shading=\"gouraud\")\n",
+ "_c = axs[0, 0].pcolormesh(x, y, np.log10(np.abs(lanczos - lag)), shading=\"gouraud\")\n",
+ "fig.colorbar(_c, ax=axs[0, :])\n",
+ "_c = axs[1, 1].pcolormesh(x, y, np.log10(rel_error), shading=\"gouraud\")\n",
+ "_c = axs[1, 0].pcolormesh(x, y, np.log10(rel_error_simple), shading=\"gouraud\")\n",
+ "fig.colorbar(_c, ax=axs[1, :])\n",
+ "_ = axs[0, 0].set_title(\"Absolute Error\")\n",
+ "_ = axs[1, 0].set_title(\"Relative Error\")\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = 0.5\n",
+ "ns = [4, 5, 5, 6, 7, 8, 8, 9, 10, 11, 11, 12] # np.arange(4, 13)\n",
+ "ms = np.arange(6, 18)\n",
+ "xi = np.logspace(0, 2, 201)[:, None]\n",
+ "lanczos = eval_lanczos([z])[0]\n",
+ "\n",
+ "_, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(12, 8))\n",
+ "ax.grid(1)\n",
+ "for n, m in zip(ns, ms):\n",
+ " zeros, weights = np.polynomial.laguerre.laggauss(n)\n",
+ " c = scipy.special.factorial(n) ** 2 / scipy.special.factorial(2 * n)\n",
+ " e = np.abs(\n",
+ " scipy.special.poch(z - 2 * n, 2 * n)\n",
+ " / scipy.special.poch(z - m, m)\n",
+ " * c\n",
+ " * xi ** (z - 2 * n + m - 1)\n",
+ " )\n",
+ " ez = np.sum(\n",
+ " scipy.special.poch(z - 2 * n, 2 * n)\n",
+ " / scipy.special.poch(z - m, m)\n",
+ " * c\n",
+ " * zeros[:, None] ** (z - 2 * n + m - 1),\n",
+ " 0,\n",
+ " )\n",
+ " lag = eval_laguerre([z], m)[0]\n",
+ " err = np.abs(lanczos - lag)\n",
+ " # print(m+z,ez)\n",
+ " # for zi,ezi in zip(z[0], ez):\n",
+ " # print(f\"{m+zi}: {ezi}\")\n",
+ " # ax.semilogy(xi, e, color=color)\n",
+ " lines = ax.loglog(xi, e, label=str(n))\n",
+ " ax.axhline(err, color=lines[0].get_color())\n",
+ " # ax.set_xticks(np.arange(xi[-1] + 1))\n",
+ " # ax.set_ylim(1e-8, 1e5)\n",
+ "_ = ax.legend()\n",
+ "# _ = ax.legend([f\"z={zi}\" for zi in z[0]])\n",
+ "# _ = [ax.axvline(x) for x in zeros]\n"
+ ]
+ }
+ ],
+ "metadata": {
+ "interpreter": {
+ "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
+ },
+ "kernelspec": {
+ "display_name": "Python 3.8.10 64-bit",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.10"
+ },
+ "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/buch/papers/laguerre/scripts/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py
new file mode 100644
index 0000000..b9088d0
--- /dev/null
+++ b/buch/papers/laguerre/scripts/laguerre_plot.py
@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""Some plots for Laguerre Polynomials."""
+
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.special as ss
+
+
+def get_ticks(start, end, step=1):
+ ticks = np.arange(start, end, step)
+ return ticks[ticks != 0]
+
+
+N = 1000
+step = 5
+t = np.linspace(-1.05, 10.5, N)[:, None]
+root = str(Path(__file__).parent)
+img_path = f"{root}/../images"
+os.makedirs(img_path, exist_ok=True)
+
+
+# fig = plt.figure(num=1, clear=True, tight_layout=True, figsize=(5.5, 3.7))
+# ax = fig.add_subplot(axes_class=AxesZero)
+fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(6, 4))
+for n in np.arange(0, 8):
+ k = np.arange(0, n + 1)[None]
+ L = np.sum((-1) ** k * ss.binom(n, k) / ss.factorial(k) * t ** k, -1)
+ ax.plot(t, L, label=f"n={n}")
+
+ax.set_xticks(get_ticks(int(t[0]), t[-1]), minor=True)
+ax.set_xticks(get_ticks(0, t[-1], step))
+ax.set_xlim(t[0], t[-1] + 0.1 * (t[1] - t[0]))
+ax.set_xlabel(r"$x$", x=1.0, labelpad=-10, rotation=0, fontsize="large")
+
+ylim = 13
+ax.set_yticks(np.arange(-ylim, ylim), minor=True)
+ax.set_yticks(np.arange(-step * (ylim // step), ylim, step))
+ax.set_ylim(-ylim, ylim)
+ax.set_ylabel(r"$y$", y=0.95, labelpad=-18, rotation=0, fontsize="large")
+
+ax.legend(ncol=2, loc=(0.125, 0.01), fontsize="large")
+
+# set the x-spine
+ax.spines[["left", "bottom"]].set_position("zero")
+ax.spines[["right", "top"]].set_visible(False)
+ax.xaxis.set_ticks_position("bottom")
+hlx = 0.4
+dx = t[-1, 0] - t[0, 0]
+dy = 2 * ylim
+hly = dy / dx * hlx
+dps = fig.dpi_scale_trans.inverted()
+bbox = ax.get_window_extent().transformed(dps)
+width, height = bbox.width, bbox.height
+
+# manual arrowhead width and length
+hw = 1.0 / 60.0 * dy
+hl = 1.0 / 30.0 * dx
+lw = 0.5 # axis line width
+ohg = 0.0 # arrow overhang
+
+# compute matching arrowhead length and width
+yhw = hw / dy * dx * height / width
+yhl = hl / dx * dy * width / height
+
+# draw x and y axis
+ax.arrow(
+ t[-1, 0] - hl,
+ 0,
+ hl,
+ 0.0,
+ fc="k",
+ ec="k",
+ lw=lw,
+ head_width=hw,
+ head_length=hl,
+ overhang=ohg,
+ length_includes_head=True,
+ clip_on=False,
+)
+
+ax.arrow(
+ 0,
+ ylim - yhl,
+ 0.0,
+ yhl,
+ fc="k",
+ ec="k",
+ lw=lw,
+ head_width=yhw,
+ head_length=yhl,
+ overhang=ohg,
+ length_includes_head=True,
+ clip_on=False,
+)
+
+fig.savefig(f"{img_path}/laguerre_polynomes.pdf")
diff --git a/buch/papers/laguerre/transformation.tex b/buch/papers/laguerre/transformation.tex
deleted file mode 100644
index 4de86b6..0000000
--- a/buch/papers/laguerre/transformation.tex
+++ /dev/null
@@ -1,31 +0,0 @@
-%
-% transformation.tex
-%
-% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
-%
-\section{Laguerre Transformation
-\label{laguerre:section:transformation}}
-\begin{align}
- L \left\{ f(x) \right\}
- =
- \tilde{f}_\alpha(n)
- =
- \int_0^\infty e^{-x} x^\alpha L_n^\alpha(x) f(x) dx
- \label{laguerre:transformation}
-\end{align}
-
-\begin{align}
- L^{-1} \left\{ \tilde{f}_\alpha(n) \right\}
- =
- f(x)
- =
- \sum_{n=0}^{\infty}
- \begin{pmatrix}
- n + \alpha \\
- n
- \end{pmatrix}^{-1}
- \frac{1}{\Gamma(\alpha + 1)}
- \tilde{f}_\alpha(n)
- L_n^\alpha(x)
- \label{laguerre:inverse_transformation}
-\end{align} \ No newline at end of file
diff --git a/buch/papers/laguerre/wasserstoff.tex b/buch/papers/laguerre/wasserstoff.tex
deleted file mode 100644
index caaa6af..0000000
--- a/buch/papers/laguerre/wasserstoff.tex
+++ /dev/null
@@ -1,29 +0,0 @@
-%
-% wasserstoff.tex
-%
-% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
-%
-\section{Radialer Schwingungsanteil eines Wasserstoffatoms
-\label{laguerre:section:radial_h_atom}}
-
-\begin{align}
- \nonumber
- - \frac{\hbar^2}{2m}
- &
- \left(
- \frac{1}{r^2} \pdv{}{r}
- \left( r^2 \pdv{}{r} \right)
- +
- \frac{1}{r^2 \sin \vartheta} \pdv{}{\vartheta}
- \left( \sin \vartheta \pdv{}{\vartheta} \right)
- +
- \frac{1}{r^2 \sin^2 \vartheta} \pdv[2]{}{\varphi}
- \right)
- u(r, \vartheta, \varphi)
- \\
- & -
- \frac{e^2}{4 \pi \epsilon_0 r} u(r, \vartheta, \varphi)
- =
- E u(r, \vartheta, \varphi)
- \label{laguerre:pdg_h_atom}
-\end{align}