aboutsummaryrefslogtreecommitdiffstats
path: root/buch/papers/laguerre
diff options
context:
space:
mode:
Diffstat (limited to 'buch/papers/laguerre')
-rw-r--r--buch/papers/laguerre/Makefile39
-rw-r--r--buch/papers/laguerre/Makefile.inc4
-rw-r--r--buch/papers/laguerre/definition.tex37
-rw-r--r--buch/papers/laguerre/eigenschaften.tex66
-rw-r--r--buch/papers/laguerre/gamma.tex518
-rw-r--r--buch/papers/laguerre/images/estimates.pdfbin0 -> 13780 bytes
-rw-r--r--buch/papers/laguerre/images/gammapaths.tex1024
-rw-r--r--buch/papers/laguerre/images/gammaplot.pdfbin0 -> 23297 bytes
-rw-r--r--buch/papers/laguerre/images/gammaplot.tex73
-rw-r--r--buch/papers/laguerre/images/integrand.pdfbin0 -> 16109 bytes
-rw-r--r--buch/papers/laguerre/images/integrand_exp.pdfbin0 -> 16951 bytes
-rw-r--r--buch/papers/laguerre/images/laguerre_poly.pdfbin0 -> 19815 bytes
-rw-r--r--buch/papers/laguerre/images/laguerre_polynomes.pdfbin16239 -> 0 bytes
-rw-r--r--buch/papers/laguerre/images/rel_error_complex.pdfbin0 -> 195590 bytes
-rw-r--r--buch/papers/laguerre/images/rel_error_mirror.pdfbin0 -> 26866 bytes
-rw-r--r--buch/papers/laguerre/images/rel_error_range.pdfbin0 -> 25105 bytes
-rw-r--r--buch/papers/laguerre/images/rel_error_shifted.pdfbin0 -> 16317 bytes
-rw-r--r--buch/papers/laguerre/images/rel_error_simple.pdfbin0 -> 23353 bytes
-rw-r--r--buch/papers/laguerre/images/targets.pdfbin0 -> 14462 bytes
-rw-r--r--buch/papers/laguerre/main.tex20
-rw-r--r--buch/papers/laguerre/packages.tex2
-rw-r--r--buch/papers/laguerre/presentation/presentation.pdfbin0 -> 394774 bytes
-rw-r--r--buch/papers/laguerre/presentation/presentation.tex134
-rw-r--r--buch/papers/laguerre/presentation/sections/gamma.tex51
-rw-r--r--buch/papers/laguerre/presentation/sections/gamma_approx.tex201
-rw-r--r--buch/papers/laguerre/presentation/sections/gaussquad.tex67
-rw-r--r--buch/papers/laguerre/presentation/sections/laguerre.tex91
-rw-r--r--buch/papers/laguerre/quadratur.tex164
-rw-r--r--buch/papers/laguerre/references.bib30
-rw-r--r--buch/papers/laguerre/scripts/estimates.py49
-rw-r--r--buch/papers/laguerre/scripts/gamma_approx.ipynb395
-rw-r--r--buch/papers/laguerre/scripts/gamma_approx.py116
-rw-r--r--buch/papers/laguerre/scripts/integrand.py42
-rw-r--r--buch/papers/laguerre/scripts/integrand_exp.py46
-rw-r--r--buch/papers/laguerre/scripts/laguerre_plot.py100
-rw-r--r--buch/papers/laguerre/scripts/laguerre_poly.py108
-rw-r--r--buch/papers/laguerre/scripts/rel_error_complex.py43
-rw-r--r--buch/papers/laguerre/scripts/rel_error_mirror.py38
-rw-r--r--buch/papers/laguerre/scripts/rel_error_range.py41
-rw-r--r--buch/papers/laguerre/scripts/rel_error_shifted.py40
-rw-r--r--buch/papers/laguerre/scripts/rel_error_simple.py41
-rw-r--r--buch/papers/laguerre/scripts/targets.py58
42 files changed, 3022 insertions, 616 deletions
diff --git a/buch/papers/laguerre/Makefile b/buch/papers/laguerre/Makefile
index 0f0985a..85a1b83 100644
--- a/buch/papers/laguerre/Makefile
+++ b/buch/papers/laguerre/Makefile
@@ -3,9 +3,42 @@
#
# (c) 2020 Prof Dr Andreas Mueller
#
+IMGFOLDER := images
+PRESFOLDER := presentation
-images: images/laguerre_polynomes.pdf
+FIGURES := \
+ images/targets.pdf \
+ images/rel_error_complex.pdf \
+ images/estimates.pdf \
+ images/integrand.pdf \
+ images/integrand_exp.pdf \
+ images/laguerre_poly.pdf \
+ images/rel_error_mirror.pdf \
+ images/rel_error_range.pdf \
+ images/rel_error_shifted.pdf \
+ images/rel_error_simple.pdf \
+ images/gammaplot.pdf
-images/laguerre_polynomes.pdf: scripts/laguerre_plot.py
- python3 scripts/laguerre_plot.py
+.PHONY: all
+all: images presentation
+.PHONY: images
+images: $(FIGURES)
+
+.PHONY: presentation
+presentation: $(PRESFOLDER)/presentation.pdf
+
+images/%.pdf images/%.pgf: scripts/%.py scripts/gamma_approx.py
+ python3 $<
+
+images/gammaplot.pdf: images/gammaplot.tex images/gammapaths.tex
+ cd $(IMGFOLDER) && latexmk -quiet -pdf gammaplot.tex
+
+$(PRESFOLDER)/%.pdf: $(PRESFOLDER)/%.tex $(FIGURES)
+ cd $(PRESFOLDER) && latexmk -quiet -pdf $(<F)
+
+.PHONY: clean
+clean:
+ rm $(FIGURES)
+ cd $(IMGFOLDER) && latexmk -C
+ cd $(PRESFOLDER) && latexmk -C \ No newline at end of file
diff --git a/buch/papers/laguerre/Makefile.inc b/buch/papers/laguerre/Makefile.inc
index 12b0935..39b5d6f 100644
--- a/buch/papers/laguerre/Makefile.inc
+++ b/buch/papers/laguerre/Makefile.inc
@@ -10,6 +10,4 @@ dependencies-laguerre = \
papers/laguerre/definition.tex \
papers/laguerre/eigenschaften.tex \
papers/laguerre/quadratur.tex \
- papers/laguerre/gamma.tex
-
-
+ papers/laguerre/gamma.tex
diff --git a/buch/papers/laguerre/definition.tex b/buch/papers/laguerre/definition.tex
index d111f6f..4729a93 100644
--- a/buch/papers/laguerre/definition.tex
+++ b/buch/papers/laguerre/definition.tex
@@ -15,14 +15,17 @@ x y''(x) + (\nu + 1 - x) y'(x) + n y(x)
n \in \mathbb{N}_0
, \quad
x \in \mathbb{R}
-.
\label{laguerre:dgl}
+.
\end{align}
+Spannenderweise wurde die verallgemeinerte Laguerre-Differentialgleichung
+zuerst von Yacovlevich Sonine (1849 - 1915) beschrieben,
+aber aufgrund ihrer Ähnlichkeit nach Laguerre benannt.
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
+weil die Lösung mit derselben Methode berechnet werden kann.
+Zusätzlich erhält man aber die Lösung für den allgmeinen Fall.
+Zur Lösung von \eqref{laguerre:dgl} verwenden wir einen
Potenzreihenansatz.
Da wir bereits wissen, dass die Lösung orthogonale Polynome sind,
erscheint dieser Ansatz sinnvoll.
@@ -44,7 +47,7 @@ y''(x)
=
\sum_{k=1}^\infty (k+1) k a_{k+1} x^{k-1}
\end{align*}
-in die Differentialgleichung ein, erhält man:
+in die Differentialgleichung ein, erhält man
\begin{align*}
\sum_{k=1}^\infty (k+1) k a_{k+1} x^k
+
@@ -118,6 +121,15 @@ 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}
+Die Laguerre-Polynome von Grad $0$ bis $7$ sind in
+Abbildung~\ref{laguerre:fig:polyeval} dargestellt.
+\begin{figure}
+\centering
+% \scalebox{0.8}{\input{papers/laguerre/images/laguerre_poly.pgf}}
+\includegraphics[width=0.9\textwidth]{papers/laguerre/images/laguerre_poly.pdf}
+\caption{Laguerre-Polynome vom Grad $0$ bis $7$}
+\label{laguerre:fig:polyeval}
+\end{figure}
\subsection{Analytische Fortsetzung}
Durch die analytische Fortsetzung erhalten wir zudem noch die zweite Lösung der
@@ -126,8 +138,10 @@ Differentialgleichung mit der Form
\Xi_n(x)
=
L_n(x) \ln(x) + \sum_{k=1}^\infty d_k x^k
+.
\end{align*}
-Nach einigen mühsamen Rechnungen,
+Nach einigen aufwändigen Rechnungen,
+% die am besten ein Computeralgebrasystem übernimmt,
die den Rahmen dieses Kapitel sprengen würden,
erhalten wir
\begin{align*}
@@ -142,16 +156,5 @@ L_n(x) \ln(x)
\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 b0cc3a3..4adbe86 100644
--- a/buch/papers/laguerre/eigenschaften.tex
+++ b/buch/papers/laguerre/eigenschaften.tex
@@ -3,35 +3,31 @@
%
% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
%
-\section{Eigenschaften
- \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.
+\section{Orthogonalität
+ \label{laguerre:section:orthogonal}}
+Im Abschnitt~\ref{laguerre:section:definition}
+haben wir die Behauptung aufgestellt,
+dass die Laguerre-Polynome orthogonal 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
+Wenn wir \eqref{laguerre:dgl} 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
+Der Beweis kann äquivalent auch über den Sturm-Liouville-Operator
\begin{align}
S
=
\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right).
\label{laguerre:slop}
\end{align}
+und den Laguerre-Operator
+\begin{align}
+\Lambda
+=
+x \frac{d}{dx^2} + (\nu + 1 -x) \frac{d}{dx}
+\end{align}
+erhalten werden,
+indem wir diese Operatoren einander gleichsetzen.
Aus der Beziehung
\begin{align}
S
@@ -49,24 +45,27 @@ Ausserdem ist ersichtlich, dass $p(x)$ die Differentialgleichung
\begin{align*}
x \frac{dp}{dx}
=
--(\nu + 1 - x) p,
+-(\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
+-\int \frac{\nu + 1 - x}{x} \, dx
+=
+-\int \frac{\nu + 1}{x} \, dx - \int 1\, dx
\\
\log p
& =
--\log \nu + 1 - x + C
+-(\nu + 1)\log x - x + c
\\
p(x)
& =
-C x^{\nu + 1} e^{-x}
+.
\end{align*}
-Eingefügt in Gleichung~\eqref{laguerre:sl-lag} erhalten wir
+Eingefügt in Gleichung~\eqref{laguerre:sl-lag} ergibt sich
\begin{align*}
\frac{C}{w(x)}
\left(
@@ -106,14 +105,9 @@ Für den rechten Rand ist die Bedingung (Gleichung~\eqref{laguerre:sllag_randb})
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}
-
+Damit können wir schlussfolgern:
+Die verallgemeinerten Laguerre-Polynome sind orthogonal
+bezüglich des Skalarproduktes auf dem Intervall $(0, \infty)$
+mit der verallgemeinerten Laguerre\--Gewichtsfunktion $w(x)=x^\nu e^{-x}$.
+Die Laguerre-Polynome ($\nu=0$) sind somit orthognal im Intervall $(0, \infty)$
+mit der Gewichtsfunktion $w(x)=e^{-x}$.
diff --git a/buch/papers/laguerre/gamma.tex b/buch/papers/laguerre/gamma.tex
index e3838b0..2e5fc06 100644
--- a/buch/papers/laguerre/gamma.tex
+++ b/buch/papers/laguerre/gamma.tex
@@ -8,8 +8,8 @@
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.
+Dabei bietet sich z.B. die Gamma-Funkion hervorragend 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
@@ -19,58 +19,514 @@ Integral der Form
\begin{align}
\Gamma(z)
& =
-\int_0^\infty t^{z-1} e^{-t} dt
+\int_0^\infty x^{z-1} e^{-x} \, dx
,
\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.
+Der Term $e^{-t}$ im Integranden und der Integrationsbereich erfüllen
+genau die Bedingungen der Laguerre-Integration.
+% Der Term $e^{-t}$ ist genau die Gewichtsfunktion der Laguerre-Integration und
+% der Definitionsbereich passt ebenfalls genau für dieses Verfahren.
+Weiter zu erwähnen ist, dass für die verallgemeinerte Laguerre-Integration die
+Gewichtsfunktion $t^\nu e^{-t}$ exakt dem Integranden für $\nu=z-1$ entspricht.
\subsubsection{Funktionalgleichung}
-Die Funktionalgleichung besagt
+Die Gamma-Funktion besitzt die gleiche Rekursionsbeziehung wie die Fakultät,
+nämlich
\begin{align}
-z \Gamma(z) = \Gamma(z+1).
+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.
+
+\subsubsection{Reflektionsformel}
+Die Reflektionsformel
+\begin{align}
+\Gamma(z) \Gamma(1 - z)
+=
+\frac{\pi}{\sin \pi z}
+,\quad
+\text{für }
+z \notin \mathbb{Z}
+\label{laguerre:gamma_refform}
+\end{align}
+stellt eine Beziehung zwischen den zwei Punkten,
+die aus der Spiegelung an der Geraden $\real z = 1/2$ hervorgehen,
+her.
+Dadurch lassen Werte der Gamma-Funktion sich für $z$ in der rechten Halbebene
+leicht in die linke Halbebene übersetzen und umgekehrt.
\subsection{Berechnung mittels Gauss-Laguerre-Quadratur}
+In den vorherigen Abschnitten haben wir gesehen,
+dass sich die Gamma-Funktion bestens für die Gauss-Laguerre-Quadratur eignet.
+Nun bieten sich uns zwei Optionen,
+diese zu berechnen:
+\begin{enumerate}
+\item Wir verwenden die verallgemeinerten Laguerre-Polynome, dann $f(x)=1$.
+\item Wir verwenden die Laguerre-Polynome, dann $f(x)=x^{z-1}$.
+\end{enumerate}
+Die erste Variante wäre optimal auf das Problem angepasst,
+allerdings müssten die Gewichte und Nullstellen für jedes $z$
+neu berechnet werden,
+da sie per Definition von $z$ abhängen.
+Dazu kommt,
+dass die Berechnung der Gewichte $A_i$ nach \cite{laguerre:Cassity1965AbcissasCA}
+\begin{align*}
+A_i
+=
+\frac{
+\Gamma(n) \Gamma(n+\nu)
+}
+{
+(n+\nu)
+\left[L_{n-1}^{\nu}(x_i)\right]^2
+}
+\end{align*}
+Evaluationen der Gamma-Funktion benötigen.
+Somit ist diese Methode eindeutig nicht geeignet für unser Vorhaben.
+
+Bei der zweiten Variante benötigen wir keine Neuberechung der Gewichte
+und Nullstellen für unterschiedliche $z$.
+In \eqref{laguerre:quadratur_gewichte} ist ersichtlich,
+dass die Gewichte einfach zu berechnen sind.
+Auch die Nullstellen können vorgängig,
+mittels eines geeigneten Verfahrens,
+aus den Polynomen bestimmt werden.
+Als problematisch könnte sich höchstens
+die zu integrierende Funktion $f(x)=x^{z-1}$ für $|z| \gg 0$ erweisen.
+Somit entscheiden wir uns aufgrund der vorherigen Punkte,
+die zweite Variante weiterzuverfolgen.
+
+\subsubsection{Direkter Ansatz}
+Wenden wir also die Gauss-Laguerre-Quadratur aus
+\eqref{laguerre:laguerrequadratur} auf die Gamma-Funktion
+\eqref{laguerre:gamma} an,
+ergibt sich
+\begin{align}
+\Gamma(z)
+\approx
+\sum_{i=1}^n x_i^{z-1} A_i.
+\label{laguerre:naive_lag}
+\end{align}
+
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/rel_error_simple.pgf}
+\includegraphics{papers/laguerre/images/rel_error_simple.pdf}
+%\vspace{-12pt}
+\caption{Relativer Fehler des direkten Ansatzes
+für verschiedene reele Werte von $z$ und Grade $n$ der Laguerre-Polynome}
+\label{laguerre:fig:rel_error_simple}
+\end{figure}
-Fehlerterm:
+Bevor wir die Gauss-Laguerre-Quadratur anwenden,
+möchten wir als ersten Schritt eine Fehlerabschätzung durchführen.
+Für den Fehlerterm \eqref{laguerre:lag_error} wird die $2n$-te Ableitung
+der zu integrierenden Funktion $f(\xi)$ benötigt.
+Für das Integral der Gamma-Funktion ergibt sich also
\begin{align*}
+\frac{d^{2n}}{d\xi^{2n}} f(\xi)
+ & =
+\frac{d^{2n}}{d\xi^{2n}} \xi^{z-1}
+\\
+ & =
+(z - 2n)_{2n} \xi^{z - 2n - 1}
+.
+\end{align*}
+Eingesetzt im Fehlerterm \eqref{laguerre:lag_error} resultiert
+\begin{align}
R_n
=
(z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z-2n-1}
-\end{align*}
+,
+\label{laguerre:gamma_err_simple}
+\end{align}
+wobei $\xi$ ein geeigneter Wert im Interval $(0, \infty)$ ist
+und $n$ der Grad des verwendeten Laguerre-Polynoms.
+Eine Fehlerabschätzung mit dem Fehlerterm stellt sich als unnütz heraus,
+da $R_n$ für $z < 2n - 1$ bei $\xi \rightarrow 0$ eine Singularität aufweist
+und für $z > 2n - 1$ bei $\xi \rightarrow \infty$ divergiert.
+Nur für den unwahrscheinlichen Fall $ z = 2n - 1$
+wäre eine Fehlerabschätzung plausibel.
-\subsubsection{Finden der optimalen Berechnungsstelle}
+Wenden wir nun also direkt die Gauss-Laguerre-Quadratur auf die Gamma-Funktion
+an.
+Dazu benötigen wir die Gewichte nach
+\eqref{laguerre:quadratur_gewichte}
+und als Stützstellen die Nullstellen des Laguerre-Polynomes $L_n$.
+Evaluieren wir den relativen Fehler unserer Approximation zeigt sich ein
+Bild wie in Abbildung~\ref{laguerre:fig:rel_error_simple}.
+Man kann sehen,
+wie der relative Fehler Nullstellen aufweist für ganzzahlige $z \leq 2n$.
+Laut der Theorie der Gauss-Quadratur auch ist das zu erwarten,
+da die Approximation via Gauss-Quadratur
+exakt ist für zu integrierende Polynome mit Grad $\leq 2n-1$
+und hinzukommt,
+dass zudem von $z$ noch $1$ abgezogen wird im Exponenten.
+Es ist ersichtlich,
+dass sich für den Polynomgrad $n$ ein Interval gibt,
+in dem der relative Fehler minimal ist.
+Links steigt der relative Fehler besonders stark an,
+während er auf der rechten Seite zu konvergieren scheint.
+Um die linke Hälfte in den Griff zu bekommen,
+könnten wir die Reflektionsformel der Gamma-Funktion ausnutzen.
+
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/rel_error_mirror.pgf}
+\includegraphics{papers/laguerre/images/rel_error_mirror.pdf}
+%\vspace{-12pt}
+\caption{Relativer Fehler des Ansatzes mit Spiegelung negativer Realwerte
+für verschiedene reele Werte von $z$ und Grade $n$ der Laguerre-Polynome}
+\label{laguerre:fig:rel_error_mirror}
+\end{figure}
+
+Spiegelt man nun $z$ mit negativem Realteil mittels der Reflektionsformel,
+ergibt sich ein stabilerer Fehler in der linken Hälfte,
+wie in Abbildung~\ref{laguerre:fig:rel_error_mirror}.
+Die Spiegelung bringt nur für wenige Werte einen,
+für praktische Anwendungen geeigneten,
+relativen Fehler.
+Wie wir aber in Abbildung~\ref{laguerre:fig:rel_error_simple} sehen konnten,
+gibt es für jeden Polynomgrad $n$ ein Intervall $[a(n), a(n) + 1]$,
+$a(n) \in \mathbb{Z}$,
+in welchem der relative Fehler minimal ist.
+Die Funktionalgleichung der Gamma-Funktion \eqref{laguerre:gamma_funktional}
+könnte uns hier helfen,
+das Problem in den Griff zu bekommen.
+
+\subsubsection{Analyse des Integranden}
+Wie wir im vorherigen Abschnitt gesehen haben,
+scheint der Integrand problematisch.
+Darum möchten wir jetzt den Integranden analysieren,
+um ihn besser verstehen zu können und
+dadurch geeignete Gegenmassnahmen zu entwickeln.
+
+% Dieser Abschnitt soll eine grafisches Verständnis dafür schaffen,
+% wieso der Integrand so problematisch ist.
+% Was das heisst sollte in Abbildung~\ref{laguerre:fig:integrand}
+% und Abbildung~\ref{laguerre:fig:integrand_exp} grafisch dargestellt werden.
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/integrand.pgf}
+\includegraphics{papers/laguerre/images/integrand.pdf}
+%\vspace{-12pt}
+\caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$}
+\label{laguerre:fig:integrand}
+\end{figure}
+
+In Abbildung~\ref{laguerre:fig:integrand} ist der Integrand $x^z$ für
+unterschiedliche Werte von $z$ dargestellt.
+Dies entspricht der zu integrierenden Funktion $f(x)$
+der Gauss-Laguerre-Quadratur für die Gamma-Funktion.
+Man erkennt,
+dass für kleine $z$ sich ein singulärer Integrand ergibt
+und auch für grosse $z$ wächst der Integrand sehr schnell an.
+Das heisst,
+die Ableitungen im Fehlerterm divergieren noch schneller
+und das wirkt sich negativ auf die Genauigkeit der Approximation aus.
+Somit lässt sich hier sagen,
+dass kleine Exponenten um $0$ genauere Resultate liefern sollten.
+
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/integrand_exp.pgf}
+\includegraphics{papers/laguerre/images/integrand_exp.pdf}
+%\vspace{-12pt}
+\caption{Integrand $x^z e^{-x}$ mit unterschiedlichen Werten für $z$}
+\label{laguerre:fig:integrand_exp}
+\end{figure}
+
+In Abbildung~\ref{laguerre:fig:integrand_exp} fügen wir
+die Dämpfung der Gewichtsfunktion $w(x)$
+der Gauss-Laguerre-Quadratur wieder hinzu
+und erhalten so wieder den kompletten Integranden $x^{z-1} e^{-x}$
+der Gamma-Funktion.
+Für negative $z$ ergeben sich immer noch Singularitäten,
+wenn $x \rightarrow 0$.
+Um $1$ wächst der Term $x^z$ schneller als die Dämpfung $e^{-x}$,
+aber für $x \rightarrow \infty$ geht der Integrand gegen $0$.
+Das führt zu glockenförmigen Kurven,
+die für grosse Exponenten $z$ nach der Stelle $x=1$ schnell anwachsen.
+Zu grosse Exponenten $z$ sind also immer noch problematisch.
+Kleine positive $z$ scheinen nun also auch zulässig zu sein.
+Damit formulieren wir die Vermutung,
+dass $a(n)$,
+welches das Intervall $[a(n), a(n) + 1]$ definiert,
+in dem der relative Fehler minimal ist,
+grösser als $0$ und kleiner als $2n-1$ ist.
+
+\subsubsection{Ansatz mit Verschiebungsterm}
+% Mittels der Funktionalgleichung \eqref{laguerre:gamma_funktional}
+% kann der Wert von $\Gamma(z)$ im Interval $z \in [a,a+1]$,
+% in dem der relative Fehler minimal ist,
+% evaluiert werden und dann mit der Funktionalgleichung zurückverschoben werden.
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
+wenn man das Problem in einem geeigneten Intervall $[a(n), a(n)+1]$,
+$a(n) \in \mathbb{Z}$,
+evaluiert und dann mit der
+Funktionalgleichung \eqref{laguerre:gamma_funktional} zurückverschiebt.
+Für dieses Vorhaben führen wir einen Verschiebungsterm $m \in \mathbb{Z}$ ein.
+Passen wir \eqref{laguerre:naive_lag}
+mit dem Verschiebungsterm $m$
+%,der $z$ and die Stelle $z_m = z + m$ verschiebt,
+an,
+ergibt sich
+\begin{align}
+\Gamma(z)
+\approx
+s(z, m) \sum_{i=1}^n x_i^{z + m - 1} A_i
+% &&
+% \text{mit }
+% s(z, m)
+% =
+% \begin{cases}
+% \displaystyle
+% \frac{1}{(z - m)_m} & \text{wenn } m \geq 0\\
+% (z + m)_{-m} & \text{wenn } m < 0
+% \end{cases}
+% .
+\label{laguerre:shifted_lag}
+\end{align}
+mit
\begin{align*}
-R_n
+s(z, m)
=
-\frac{(z - 2n)_{2n}}{(z - m)_m} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1}
+\begin{cases}
+\displaystyle
+\frac{1}{(z)_m} & \text{wenn } m \geq 0 \\
+(z + m)_{-m} & \text{wenn } m < 0
+\end{cases}
.
\end{align*}
-{
-\large \color{red}
-TODO:
-Geeignete Minimierung für Fehler finden, so dass sie mit den emprisich
-bestimmen optimalen Punkten übereinstimmen.
-}
+\subsubsection{Finden der optimalen Berechnungsstelle}
+Um die optimale Stelle $z^*(n) \in \left[a(n), a(n) + 1\right]$,
+$z^*(n) \in \mathbb{R}$,
+zu finden,
+erweitern wir denn Fehlerterm \eqref{laguerre:gamma_err_simple}
+und erhalten
+\begin{align}
+R_{n,m}(\xi)
+=
+s(z, m) \cdot (z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1}
+,\quad
+\text{für }
+\xi \in (0, \infty)
+\label{laguerre:gamma_err_shifted}
+.
+\end{align}
+
+\begin{figure}
+\centering
+\includegraphics{papers/laguerre/images/targets.pdf}
+% %\vspace{-12pt}
+\caption{$a$ in Abhängigkeit von $z$ und $n$}
+\label{laguerre:fig:targets}
+\end{figure}
+% wobei ist
+% mit $z^*(n) \in \mathbb{R}$ wollen wir finden,
+% in dem wir den Fehlerterm \eqref{laguerre:lag_error} anpassen
+% und in einem nächsten Schritt minimieren.
+% Zudem nehmen wir an,
+% dass $z < z^*(n)$ ist.
+% Wir fügen einen Verschiebungsterm um $m \in \mathbb{N}$ Stellen ein,
+% daraus folgt
+%
+% Damit wir den idealen Verschiebungsterm $m^*$ finden können,
+% müssen wir mittels des Fehlerterms \eqref{laguerre:gamma_err_shifted}
+% ein Optimierungsproblem
+%
+% Das Optimierungsproblem daraus lässt sich als
+Daraus formulieren wir das Optimierungproblem
+\begin{align*}
+m^*
+=
+\operatorname*{argmin}_m \max_\xi R_{n,m}(\xi)
+.
+\end{align*}
+Allerdings ist die Funktion $R_{n,m}(\xi)$ unbeschränkt und
+hat die gleichen Probleme wie die Fehlerabschätzung des direkten Ansatzes.
+Dazu müssten wir $\xi$ versuchen,
+unter Kontrolle zu bringen,
+was ein äussersts schwieriges Unterfangen zu sein scheint.
+Da die Gauss-Quadratur aber sowieso
+nur wirklich praktisch sinnvoll für kleine $n$ ist,
+können die Intervalle $[a(n), a(n)+1]$ empirisch gesucht werden.
+
+Wir bestimmen nun die optimalen Verschiebungsterme empirisch
+für $n = 2,\ldots, 12$ im Intervall $z \in (0, 1)$,
+da $z$ sowieso um den Term $m$ verschoben wird,
+reicht die $m^*$ nur in diesem Intervall zu analysieren.
+In Abbildung~\ref{laguerre:fig:targets} sind die empirisch bestimmten $m^*$
+abhängig von $z$ und $n$ dargestellt.
+In $n$-Richtung lässt sich eine klare lineare Abhängigkeit erkennen und
+die Beziehung zu $z$ ist negativ,
+d.h. wenn $z$ grösser ist, wird $m^*$ kleiner.
+Allerdings ist die genaue Beziehung zu $z$
+aus dieser Grafik nicht offensichtlich,
+aber sie scheint regelmässig zu sein.
+Es lässt die Vermutung aufkommen,
+dass die Restriktion von $m^* \in \mathbb{Z}$ Rundungsprobleme verursacht.
+Wir versuchen,
+dieses Problem via lineare Regression und geeignete Rundung zu beheben.
+Den linearen Regressor
+\begin{align*}
+\hat{m}
+=
+\alpha n + \beta
+\end{align*}
+machen wir nur abhängig von $n$
+in dem wir den Mittelwert $\overline{m}$ von $m^*$ über $z$ berechnen.
+
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/estimates.pgf}
+\includegraphics{papers/laguerre/images/estimates.pdf}
+%\vspace{-12pt}
+\caption{Schätzung Mittelwert von $m$ und Fehler}
+\label{laguerre:fig:schaetzung}
+\end{figure}
-\subsection{Resultate}
+In Abbildung~\ref{laguerre:fig:schaetzung} sind die Resultate
+der linearen Regression aufgezeigt mit $\alpha = 1.34094$ und $\beta =
+0.854093$.
+Die lineare Beziehung ist ganz klar ersichtlich und der Fit scheint zu genügen.
+Der optimale Verschiebungsterm kann nun mit
+\begin{align*}
+m^*
+\approx
+\lceil \hat{m} - z \rceil
+=
+\lceil \alpha n + \beta - z \rceil
+\end{align*}
+% kann nun mit dem linearen Regressor und $z$
+gefunden werden.
+
+\subsubsection{Evaluation des Schätzers}
+In einem ersten Schritt möchten wir analysieren,
+wie gut die Abschätzung des optimalen Verschiebungsterms ist.
+Dazu bestimmen wir den relativen Fehler für verschiedene Verschiebungsterme $m$
+rund um $m^*$ bei gegebenem Polynomgrad $n = 8$ für $z \in (0, 1)$.
+Abbildung~\ref{laguerre:fig:rel_error_shifted} sind die relativen Fehler
+der Approximation dargestellt.
+Man kann deutlich sehen,
+dass der relative Fehler anwächst,
+je weiter der Verschiebungsterm vom idealen Wert abweicht.
+Zudem scheint der Schätzer den optimalen Verschiebungsterm gut zu bestimmen,
+da der Schätzer zuerst der grünen Linie folgt und
+dann beim Übergang auf die orange Linie wechselt.
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/rel_error_shifted.pgf}
+\includegraphics{papers/laguerre/images/rel_error_shifted.pdf}
+%\vspace{-12pt}
+\caption{Relativer Fehler des Ansatzes mit Verschiebungsterm
+für verschiedene reele Werte von $z$ und Verschiebungsterme $m$.
+Das verwendete Laguerre-Polynom besitzt den Grad $n = 8$.
+$m^*$ bezeichnet hier den optimalen Verschiebungsterm.}
+\label{laguerre:fig:rel_error_shifted}
+\end{figure}
+
+\subsubsection{Resultate}
+Das Verfahren scheint für den Grad $n=8$ und $z \in (0,1)$ gut zu funktioneren.
+Es stellt sich nun die Frage,
+wie der relative Fehler sich für verschiedene $z$ und $n$ verhält.
+In Abbildung~\ref{laguerre:fig:rel_error_range} sind die relativen Fehler für
+unterschiedliche $n$ dargestellt.
+Der relative Fehler scheint immer noch Nullstellen aufzuweisen
+für ganzzahlige $z$.
+Durch das Verschieben ergibt sich jetzt aber,
+wie zu erwarten war,
+ein periodischer relativer Fehler mit einer Periodendauer von $1$.
+Zudem lässt sich erkennen,
+dass der Fehler abhängig von der Ordnung $n$
+des verwendeten Laguerre-Polynoms ist.
+Wenn der Grad $n$ um $1$ erhöht wird,
+verbessert sich die Genauigkeit des Resultats um etwa eine signifikante Stelle.
+
+In Abbildung~\ref{laguerre:fig:rel_error_complex}
+ist der Betrag des relativen Fehlers in der komplexen Ebene dargestellt.
+Je stärker der Imaginäranteil von $z$ von $0$ abweicht,
+umso schlechter wird die Genauigkeit der Approximation.
+Das erstaunt nicht weiter,
+da die Gauss-Quadratur eigentlich nur für reelle Zahlen definiert ist.
+Wenn der Imaginäranteil von $z$ ungefähr $0$ ist,
+lässt sich das gleiche Bild beobachten wie in
+Abbildung~\ref{laguerre:fig:rel_error_range}.
+
+\begin{figure}
+\centering
+% \input{papers/laguerre/images/rel_error_range.pgf}
+\includegraphics{papers/laguerre/images/rel_error_range.pdf}
+%\vspace{-12pt}
+\caption{Relativer Fehler des Ansatzes mit optimalen Verschiebungsterm
+für verschiedene reele Werte von $z$ und Laguerre-Polynome vom Grad $n$}
+\label{laguerre:fig:rel_error_range}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics{papers/laguerre/images/rel_error_complex.pdf}
+%\vspace{-12pt}
+\caption{Absolutwert des relativen Fehlers in der komplexen Ebene}
+\label{laguerre:fig:rel_error_complex}
+\end{figure}
+
+\subsubsection{Vergleich mit Lanczos-Methode}
+Nun stellt sich die Frage,
+wie das in diesem Abschnitt beschriebene Approximationsverfahren
+der Gamma-Funktion sich gegenüber den üblichen Approximationsverfahren schlägt.
+Eine häufig verwendete Methode ist die Lanczos-Approximation,
+welche gegeben ist durch
+\begin{align}
+\Gamma(z + 1)
+\approx
+\sqrt{2\pi} \left( z + \sigma + \frac{1}{2} \right)^{z + 1/2}
+e^{-(z + \sigma + 1/2)} \sum_{k=0}^n g_k H_k(z)
+,
+\end{align}
+wobei
+\begin{align*}
+g_k = \frac{e^\sigma \varepsilon_k (-1)^k}{\sqrt{2\pi}}
+\sum_{r=0}^k (-1)^r \, \binom{k}{r} \, (k)_r
+\left( \frac{e}{r + \sigma + \frac{1}{2}}\right)^{r + 1/2}
+,
+\end{align*}
+\begin{align*}
+\varepsilon_k
+=
+\begin{cases}
+1 & \text{für } k = 0 \\
+2 & \text{sonst}
+\end{cases}
+\quad \text{und}\quad
+H_k(z)
+=
+\frac{(-1)^k (-z)_k}{(z+1)_k}
+\end{align*}
+mit $H_0 = 1$ und $\sum_0^n g_k = 1$ (siehe \cite{laguerre:lanczos}).
+Diese Methode wurde zum Beispiel in
+{\em GNU Scientific Library}, {\em Boost}, {\em CPython} und
+{\em musl} implementiert.
+Diese Methode erreicht für $n = 7$ typischerweise Genauigkeit von $13$
+korrekten, signifikanten Stellen für reele Argumente.
+Zum Vergleich: die vorgestellte Methode erreicht für $n = 7$
+eine minimale Genauigkeit von $6$ korrekten, signifikanten Stellen
+für reele Argumente.
+Das Resultat ist etwas enttäuschend,
+aber nicht unerwartet,
+da die Lanczos-Methode spezifisch auf dieses Problem zugeschnitten ist und
+unsere Methode eine erweiterte allgemeine Methode ist.
+Was die Komplexität der Berechnungen im Betrieb angeht,
+ist die Gauss-Laguerre-Quadratur wesentlich ressourcensparender,
+weil sie nur aus $n$ Funktionsevaluationen,
+wenigen Multiplikationen und Additionen besteht.
+Demzufolge könnte diese Methode Anwendung in Systemen mit wenig Rechenleistung
+und/oder knappen Energieressourcen finden. \ No newline at end of file
diff --git a/buch/papers/laguerre/images/estimates.pdf b/buch/papers/laguerre/images/estimates.pdf
new file mode 100644
index 0000000..bd995de
--- /dev/null
+++ b/buch/papers/laguerre/images/estimates.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/gammapaths.tex b/buch/papers/laguerre/images/gammapaths.tex
new file mode 100644
index 0000000..efa0863
--- /dev/null
+++ b/buch/papers/laguerre/images/gammapaths.tex
@@ -0,0 +1,1024 @@
+\def\gammaplus{({\dx*0.0190},{\dy*52.0728})
+ -- ({\dx*0.0200},{\dy*49.4422})
+ -- ({\dx*0.0400},{\dy*24.4610})
+ -- ({\dx*0.0600},{\dy*16.1457})
+ -- ({\dx*0.0800},{\dy*11.9966})
+ -- ({\dx*0.1000},{\dy*9.5135})
+ -- ({\dx*0.1200},{\dy*7.8633})
+ -- ({\dx*0.1400},{\dy*6.6887})
+ -- ({\dx*0.1600},{\dy*5.8113})
+ -- ({\dx*0.1800},{\dy*5.1318})
+ -- ({\dx*0.2000},{\dy*4.5908})
+ -- ({\dx*0.2200},{\dy*4.1505})
+ -- ({\dx*0.2400},{\dy*3.7855})
+ -- ({\dx*0.2600},{\dy*3.4785})
+ -- ({\dx*0.2800},{\dy*3.2169})
+ -- ({\dx*0.3000},{\dy*2.9916})
+ -- ({\dx*0.3200},{\dy*2.7958})
+ -- ({\dx*0.3400},{\dy*2.6242})
+ -- ({\dx*0.3600},{\dy*2.4727})
+ -- ({\dx*0.3800},{\dy*2.3383})
+ -- ({\dx*0.4000},{\dy*2.2182})
+ -- ({\dx*0.4200},{\dy*2.1104})
+ -- ({\dx*0.4400},{\dy*2.0132})
+ -- ({\dx*0.4600},{\dy*1.9252})
+ -- ({\dx*0.4800},{\dy*1.8453})
+ -- ({\dx*0.5000},{\dy*1.7725})
+ -- ({\dx*0.5200},{\dy*1.7058})
+ -- ({\dx*0.5400},{\dy*1.6448})
+ -- ({\dx*0.5600},{\dy*1.5886})
+ -- ({\dx*0.5800},{\dy*1.5369})
+ -- ({\dx*0.6000},{\dy*1.4892})
+ -- ({\dx*0.6200},{\dy*1.4450})
+ -- ({\dx*0.6400},{\dy*1.4041})
+ -- ({\dx*0.6600},{\dy*1.3662})
+ -- ({\dx*0.6800},{\dy*1.3309})
+ -- ({\dx*0.7000},{\dy*1.2981})
+ -- ({\dx*0.7200},{\dy*1.2675})
+ -- ({\dx*0.7400},{\dy*1.2390})
+ -- ({\dx*0.7600},{\dy*1.2123})
+ -- ({\dx*0.7800},{\dy*1.1875})
+ -- ({\dx*0.8000},{\dy*1.1642})
+ -- ({\dx*0.8200},{\dy*1.1425})
+ -- ({\dx*0.8400},{\dy*1.1222})
+ -- ({\dx*0.8600},{\dy*1.1031})
+ -- ({\dx*0.8800},{\dy*1.0853})
+ -- ({\dx*0.9000},{\dy*1.0686})
+ -- ({\dx*0.9200},{\dy*1.0530})
+ -- ({\dx*0.9400},{\dy*1.0384})
+ -- ({\dx*0.9600},{\dy*1.0247})
+ -- ({\dx*0.9800},{\dy*1.0119})
+ -- ({\dx*1.0000},{\dy*1.0000})
+ -- ({\dx*1.0200},{\dy*0.9888})
+ -- ({\dx*1.0400},{\dy*0.9784})
+ -- ({\dx*1.0600},{\dy*0.9687})
+ -- ({\dx*1.0800},{\dy*0.9597})
+ -- ({\dx*1.1000},{\dy*0.9514})
+ -- ({\dx*1.1200},{\dy*0.9436})
+ -- ({\dx*1.1400},{\dy*0.9364})
+ -- ({\dx*1.1600},{\dy*0.9298})
+ -- ({\dx*1.1800},{\dy*0.9237})
+ -- ({\dx*1.2000},{\dy*0.9182})
+ -- ({\dx*1.2200},{\dy*0.9131})
+ -- ({\dx*1.2400},{\dy*0.9085})
+ -- ({\dx*1.2600},{\dy*0.9044})
+ -- ({\dx*1.2800},{\dy*0.9007})
+ -- ({\dx*1.3000},{\dy*0.8975})
+ -- ({\dx*1.3200},{\dy*0.8946})
+ -- ({\dx*1.3400},{\dy*0.8922})
+ -- ({\dx*1.3600},{\dy*0.8902})
+ -- ({\dx*1.3800},{\dy*0.8885})
+ -- ({\dx*1.4000},{\dy*0.8873})
+ -- ({\dx*1.4200},{\dy*0.8864})
+ -- ({\dx*1.4400},{\dy*0.8858})
+ -- ({\dx*1.4600},{\dy*0.8856})
+ -- ({\dx*1.4800},{\dy*0.8857})
+ -- ({\dx*1.5000},{\dy*0.8862})
+ -- ({\dx*1.5200},{\dy*0.8870})
+ -- ({\dx*1.5400},{\dy*0.8882})
+ -- ({\dx*1.5600},{\dy*0.8896})
+ -- ({\dx*1.5800},{\dy*0.8914})
+ -- ({\dx*1.6000},{\dy*0.8935})
+ -- ({\dx*1.6200},{\dy*0.8959})
+ -- ({\dx*1.6400},{\dy*0.8986})
+ -- ({\dx*1.6600},{\dy*0.9017})
+ -- ({\dx*1.6800},{\dy*0.9050})
+ -- ({\dx*1.7000},{\dy*0.9086})
+ -- ({\dx*1.7200},{\dy*0.9126})
+ -- ({\dx*1.7400},{\dy*0.9168})
+ -- ({\dx*1.7600},{\dy*0.9214})
+ -- ({\dx*1.7800},{\dy*0.9262})
+ -- ({\dx*1.8000},{\dy*0.9314})
+ -- ({\dx*1.8200},{\dy*0.9368})
+ -- ({\dx*1.8400},{\dy*0.9426})
+ -- ({\dx*1.8600},{\dy*0.9487})
+ -- ({\dx*1.8800},{\dy*0.9551})
+ -- ({\dx*1.9000},{\dy*0.9618})
+ -- ({\dx*1.9200},{\dy*0.9688})
+ -- ({\dx*1.9400},{\dy*0.9761})
+ -- ({\dx*1.9600},{\dy*0.9837})
+ -- ({\dx*1.9800},{\dy*0.9917})
+ -- ({\dx*2.0000},{\dy*1.0000})
+ -- ({\dx*2.0200},{\dy*1.0086})
+ -- ({\dx*2.0400},{\dy*1.0176})
+ -- ({\dx*2.0600},{\dy*1.0269})
+ -- ({\dx*2.0800},{\dy*1.0365})
+ -- ({\dx*2.1000},{\dy*1.0465})
+ -- ({\dx*2.1200},{\dy*1.0568})
+ -- ({\dx*2.1400},{\dy*1.0675})
+ -- ({\dx*2.1600},{\dy*1.0786})
+ -- ({\dx*2.1800},{\dy*1.0900})
+ -- ({\dx*2.2000},{\dy*1.1018})
+ -- ({\dx*2.2200},{\dy*1.1140})
+ -- ({\dx*2.2400},{\dy*1.1266})
+ -- ({\dx*2.2600},{\dy*1.1395})
+ -- ({\dx*2.2800},{\dy*1.1529})
+ -- ({\dx*2.3000},{\dy*1.1667})
+ -- ({\dx*2.3200},{\dy*1.1809})
+ -- ({\dx*2.3400},{\dy*1.1956})
+ -- ({\dx*2.3600},{\dy*1.2107})
+ -- ({\dx*2.3800},{\dy*1.2262})
+ -- ({\dx*2.4000},{\dy*1.2422})
+ -- ({\dx*2.4200},{\dy*1.2586})
+ -- ({\dx*2.4400},{\dy*1.2756})
+ -- ({\dx*2.4600},{\dy*1.2930})
+ -- ({\dx*2.4800},{\dy*1.3109})
+ -- ({\dx*2.5000},{\dy*1.3293})
+ -- ({\dx*2.5200},{\dy*1.3483})
+ -- ({\dx*2.5400},{\dy*1.3678})
+ -- ({\dx*2.5600},{\dy*1.3878})
+ -- ({\dx*2.5800},{\dy*1.4084})
+ -- ({\dx*2.6000},{\dy*1.4296})
+ -- ({\dx*2.6200},{\dy*1.4514})
+ -- ({\dx*2.6400},{\dy*1.4738})
+ -- ({\dx*2.6600},{\dy*1.4968})
+ -- ({\dx*2.6800},{\dy*1.5204})
+ -- ({\dx*2.7000},{\dy*1.5447})
+ -- ({\dx*2.7200},{\dy*1.5696})
+ -- ({\dx*2.7400},{\dy*1.5953})
+ -- ({\dx*2.7600},{\dy*1.6216})
+ -- ({\dx*2.7800},{\dy*1.6487})
+ -- ({\dx*2.8000},{\dy*1.6765})
+ -- ({\dx*2.8200},{\dy*1.7051})
+ -- ({\dx*2.8400},{\dy*1.7344})
+ -- ({\dx*2.8600},{\dy*1.7646})
+ -- ({\dx*2.8800},{\dy*1.7955})
+ -- ({\dx*2.9000},{\dy*1.8274})
+ -- ({\dx*2.9200},{\dy*1.8600})
+ -- ({\dx*2.9400},{\dy*1.8936})
+ -- ({\dx*2.9600},{\dy*1.9281})
+ -- ({\dx*2.9800},{\dy*1.9636})
+ -- ({\dx*3.0000},{\dy*2.0000})
+ -- ({\dx*3.0200},{\dy*2.0374})
+ -- ({\dx*3.0400},{\dy*2.0759})
+ -- ({\dx*3.0600},{\dy*2.1153})
+ -- ({\dx*3.0800},{\dy*2.1559})
+ -- ({\dx*3.1000},{\dy*2.1976})
+ -- ({\dx*3.1200},{\dy*2.2405})
+ -- ({\dx*3.1400},{\dy*2.2845})
+ -- ({\dx*3.1600},{\dy*2.3297})
+ -- ({\dx*3.1800},{\dy*2.3762})
+ -- ({\dx*3.2000},{\dy*2.4240})
+ -- ({\dx*3.2200},{\dy*2.4731})
+ -- ({\dx*3.2400},{\dy*2.5235})
+ -- ({\dx*3.2600},{\dy*2.5754})
+ -- ({\dx*3.2800},{\dy*2.6287})
+ -- ({\dx*3.3000},{\dy*2.6834})
+ -- ({\dx*3.3200},{\dy*2.7397})
+ -- ({\dx*3.3400},{\dy*2.7976})
+ -- ({\dx*3.3600},{\dy*2.8571})
+ -- ({\dx*3.3800},{\dy*2.9183})
+ -- ({\dx*3.4000},{\dy*2.9812})
+ -- ({\dx*3.4200},{\dy*3.0459})
+ -- ({\dx*3.4400},{\dy*3.1124})
+ -- ({\dx*3.4600},{\dy*3.1807})
+ -- ({\dx*3.4800},{\dy*3.2510})
+ -- ({\dx*3.5000},{\dy*3.3234})
+ -- ({\dx*3.5200},{\dy*3.3977})
+ -- ({\dx*3.5400},{\dy*3.4742})
+ -- ({\dx*3.5600},{\dy*3.5529})
+ -- ({\dx*3.5800},{\dy*3.6338})
+ -- ({\dx*3.6000},{\dy*3.7170})
+ -- ({\dx*3.6200},{\dy*3.8027})
+ -- ({\dx*3.6400},{\dy*3.8908})
+ -- ({\dx*3.6600},{\dy*3.9814})
+ -- ({\dx*3.6800},{\dy*4.0747})
+ -- ({\dx*3.7000},{\dy*4.1707})
+ -- ({\dx*3.7200},{\dy*4.2694})
+ -- ({\dx*3.7400},{\dy*4.3711})
+ -- ({\dx*3.7600},{\dy*4.4757})
+ -- ({\dx*3.7800},{\dy*4.5833})
+ -- ({\dx*3.8000},{\dy*4.6942})
+ -- ({\dx*3.8200},{\dy*4.8083})
+ -- ({\dx*3.8400},{\dy*4.9257})
+ -- ({\dx*3.8600},{\dy*5.0466})
+ -- ({\dx*3.8800},{\dy*5.1711})
+ -- ({\dx*3.9000},{\dy*5.2993})
+ -- ({\dx*3.9200},{\dy*5.4313})
+ -- ({\dx*3.9400},{\dy*5.5673})
+ -- ({\dx*3.9600},{\dy*5.7073})
+ -- ({\dx*3.9800},{\dy*5.8515})
+ -- ({\dx*4.0000},{\dy*6.0000})
+ -- ({\dx*4.0200},{\dy*6.1530})
+ -- ({\dx*4.0400},{\dy*6.3106})
+ -- ({\dx*4.0600},{\dy*6.4730})
+ -- ({\dx*4.0800},{\dy*6.6403})
+ -- ({\dx*4.0810},{\dy*6.6488})}
+\def\gammaone{({\dx*-0.9810},{\dy*-53.0814})
+ -- ({\dx*-0.9800},{\dy*-50.4512})
+ -- ({\dx*-0.9600},{\dy*-25.4802})
+ -- ({\dx*-0.9400},{\dy*-17.1763})
+ -- ({\dx*-0.9200},{\dy*-13.0397})
+ -- ({\dx*-0.9000},{\dy*-10.5706})
+ -- ({\dx*-0.8800},{\dy*-8.9355})
+ -- ({\dx*-0.8600},{\dy*-7.7775})
+ -- ({\dx*-0.8400},{\dy*-6.9182})
+ -- ({\dx*-0.8200},{\dy*-6.2583})
+ -- ({\dx*-0.8000},{\dy*-5.7386})
+ -- ({\dx*-0.7800},{\dy*-5.3211})
+ -- ({\dx*-0.7600},{\dy*-4.9809})
+ -- ({\dx*-0.7400},{\dy*-4.7006})
+ -- ({\dx*-0.7200},{\dy*-4.4678})
+ -- ({\dx*-0.7000},{\dy*-4.2737})
+ -- ({\dx*-0.6800},{\dy*-4.1114})
+ -- ({\dx*-0.6600},{\dy*-3.9760})
+ -- ({\dx*-0.6400},{\dy*-3.8636})
+ -- ({\dx*-0.6200},{\dy*-3.7714})
+ -- ({\dx*-0.6000},{\dy*-3.6969})
+ -- ({\dx*-0.5800},{\dy*-3.6386})
+ -- ({\dx*-0.5600},{\dy*-3.5950})
+ -- ({\dx*-0.5400},{\dy*-3.5652})
+ -- ({\dx*-0.5200},{\dy*-3.5487})
+ -- ({\dx*-0.5000},{\dy*-3.5449})
+ -- ({\dx*-0.4800},{\dy*-3.5538})
+ -- ({\dx*-0.4600},{\dy*-3.5756})
+ -- ({\dx*-0.4400},{\dy*-3.6105})
+ -- ({\dx*-0.4200},{\dy*-3.6594})
+ -- ({\dx*-0.4000},{\dy*-3.7230})
+ -- ({\dx*-0.3800},{\dy*-3.8027})
+ -- ({\dx*-0.3600},{\dy*-3.9004})
+ -- ({\dx*-0.3400},{\dy*-4.0181})
+ -- ({\dx*-0.3200},{\dy*-4.1590})
+ -- ({\dx*-0.3000},{\dy*-4.3269})
+ -- ({\dx*-0.2800},{\dy*-4.5267})
+ -- ({\dx*-0.2600},{\dy*-4.7652})
+ -- ({\dx*-0.2400},{\dy*-5.0514})
+ -- ({\dx*-0.2200},{\dy*-5.3976})
+ -- ({\dx*-0.2000},{\dy*-5.8211})
+ -- ({\dx*-0.1800},{\dy*-6.3472})
+ -- ({\dx*-0.1600},{\dy*-7.0135})
+ -- ({\dx*-0.1400},{\dy*-7.8795})
+ -- ({\dx*-0.1200},{\dy*-9.0442})
+ -- ({\dx*-0.1000},{\dy*-10.6863})
+ -- ({\dx*-0.0800},{\dy*-13.1627})
+ -- ({\dx*-0.0600},{\dy*-17.3067})
+ -- ({\dx*-0.0400},{\dy*-25.6183})
+ -- ({\dx*-0.0200},{\dy*-50.5974})
+ -- ({\dx*-0.0190},{\dy*-53.2279})}
+\def\gammatwo{({\dx*-1.9810},{\dy*26.7952})
+ -- ({\dx*-1.9800},{\dy*25.4804})
+ -- ({\dx*-1.9600},{\dy*13.0001})
+ -- ({\dx*-1.9400},{\dy*8.8538})
+ -- ({\dx*-1.9200},{\dy*6.7915})
+ -- ({\dx*-1.9000},{\dy*5.5635})
+ -- ({\dx*-1.8800},{\dy*4.7529})
+ -- ({\dx*-1.8600},{\dy*4.1815})
+ -- ({\dx*-1.8400},{\dy*3.7599})
+ -- ({\dx*-1.8200},{\dy*3.4386})
+ -- ({\dx*-1.8000},{\dy*3.1881})
+ -- ({\dx*-1.7800},{\dy*2.9894})
+ -- ({\dx*-1.7600},{\dy*2.8301})
+ -- ({\dx*-1.7400},{\dy*2.7015})
+ -- ({\dx*-1.7200},{\dy*2.5976})
+ -- ({\dx*-1.7000},{\dy*2.5139})
+ -- ({\dx*-1.6800},{\dy*2.4473})
+ -- ({\dx*-1.6600},{\dy*2.3952})
+ -- ({\dx*-1.6400},{\dy*2.3559})
+ -- ({\dx*-1.6200},{\dy*2.3280})
+ -- ({\dx*-1.6000},{\dy*2.3106})
+ -- ({\dx*-1.5800},{\dy*2.3029})
+ -- ({\dx*-1.5600},{\dy*2.3045})
+ -- ({\dx*-1.5400},{\dy*2.3151})
+ -- ({\dx*-1.5200},{\dy*2.3346})
+ -- ({\dx*-1.5000},{\dy*2.3633})
+ -- ({\dx*-1.4800},{\dy*2.4012})
+ -- ({\dx*-1.4600},{\dy*2.4490})
+ -- ({\dx*-1.4400},{\dy*2.5073})
+ -- ({\dx*-1.4200},{\dy*2.5770})
+ -- ({\dx*-1.4000},{\dy*2.6593})
+ -- ({\dx*-1.3800},{\dy*2.7556})
+ -- ({\dx*-1.3600},{\dy*2.8679})
+ -- ({\dx*-1.3400},{\dy*2.9986})
+ -- ({\dx*-1.3200},{\dy*3.1508})
+ -- ({\dx*-1.3000},{\dy*3.3283})
+ -- ({\dx*-1.2800},{\dy*3.5365})
+ -- ({\dx*-1.2600},{\dy*3.7819})
+ -- ({\dx*-1.2400},{\dy*4.0737})
+ -- ({\dx*-1.2200},{\dy*4.4243})
+ -- ({\dx*-1.2000},{\dy*4.8510})
+ -- ({\dx*-1.1800},{\dy*5.3790})
+ -- ({\dx*-1.1600},{\dy*6.0461})
+ -- ({\dx*-1.1400},{\dy*6.9118})
+ -- ({\dx*-1.1200},{\dy*8.0752})
+ -- ({\dx*-1.1000},{\dy*9.7148})
+ -- ({\dx*-1.0800},{\dy*12.1877})
+ -- ({\dx*-1.0600},{\dy*16.3271})
+ -- ({\dx*-1.0400},{\dy*24.6330})
+ -- ({\dx*-1.0200},{\dy*49.6053})
+ -- ({\dx*-1.0190},{\dy*52.2354})}
+\def\gammathree{({\dx*-2.9810},{\dy*-8.9887})
+ -- ({\dx*-2.9800},{\dy*-8.5505})
+ -- ({\dx*-2.9600},{\dy*-4.3919})
+ -- ({\dx*-2.9400},{\dy*-3.0115})
+ -- ({\dx*-2.9200},{\dy*-2.3259})
+ -- ({\dx*-2.9000},{\dy*-1.9184})
+ -- ({\dx*-2.8800},{\dy*-1.6503})
+ -- ({\dx*-2.8600},{\dy*-1.4621})
+ -- ({\dx*-2.8400},{\dy*-1.3239})
+ -- ({\dx*-2.8200},{\dy*-1.2194})
+ -- ({\dx*-2.8000},{\dy*-1.1386})
+ -- ({\dx*-2.7800},{\dy*-1.0753})
+ -- ({\dx*-2.7600},{\dy*-1.0254})
+ -- ({\dx*-2.7400},{\dy*-0.9859})
+ -- ({\dx*-2.7200},{\dy*-0.9550})
+ -- ({\dx*-2.7000},{\dy*-0.9311})
+ -- ({\dx*-2.6800},{\dy*-0.9132})
+ -- ({\dx*-2.6600},{\dy*-0.9004})
+ -- ({\dx*-2.6400},{\dy*-0.8924})
+ -- ({\dx*-2.6200},{\dy*-0.8886})
+ -- ({\dx*-2.6000},{\dy*-0.8887})
+ -- ({\dx*-2.5800},{\dy*-0.8926})
+ -- ({\dx*-2.5600},{\dy*-0.9002})
+ -- ({\dx*-2.5400},{\dy*-0.9115})
+ -- ({\dx*-2.5200},{\dy*-0.9264})
+ -- ({\dx*-2.5000},{\dy*-0.9453})
+ -- ({\dx*-2.4800},{\dy*-0.9682})
+ -- ({\dx*-2.4600},{\dy*-0.9955})
+ -- ({\dx*-2.4400},{\dy*-1.0276})
+ -- ({\dx*-2.4200},{\dy*-1.0649})
+ -- ({\dx*-2.4000},{\dy*-1.1080})
+ -- ({\dx*-2.3800},{\dy*-1.1578})
+ -- ({\dx*-2.3600},{\dy*-1.2152})
+ -- ({\dx*-2.3400},{\dy*-1.2815})
+ -- ({\dx*-2.3200},{\dy*-1.3581})
+ -- ({\dx*-2.3000},{\dy*-1.4471})
+ -- ({\dx*-2.2800},{\dy*-1.5511})
+ -- ({\dx*-2.2600},{\dy*-1.6734})
+ -- ({\dx*-2.2400},{\dy*-1.8186})
+ -- ({\dx*-2.2200},{\dy*-1.9929})
+ -- ({\dx*-2.2000},{\dy*-2.2050})
+ -- ({\dx*-2.1800},{\dy*-2.4674})
+ -- ({\dx*-2.1600},{\dy*-2.7991})
+ -- ({\dx*-2.1400},{\dy*-3.2298})
+ -- ({\dx*-2.1200},{\dy*-3.8091})
+ -- ({\dx*-2.1000},{\dy*-4.6261})
+ -- ({\dx*-2.0800},{\dy*-5.8595})
+ -- ({\dx*-2.0600},{\dy*-7.9258})
+ -- ({\dx*-2.0400},{\dy*-12.0750})
+ -- ({\dx*-2.0200},{\dy*-24.5571})
+ -- ({\dx*-2.0190},{\dy*-25.8719})}
+\def\gammafour{({\dx*-3.9950},{\dy*8.3966})
+ -- ({\dx*-3.9800},{\dy*2.1484})
+ -- ({\dx*-3.9600},{\dy*1.1091})
+ -- ({\dx*-3.9400},{\dy*0.7643})
+ -- ({\dx*-3.9200},{\dy*0.5933})
+ -- ({\dx*-3.9000},{\dy*0.4919})
+ -- ({\dx*-3.8800},{\dy*0.4253})
+ -- ({\dx*-3.8600},{\dy*0.3788})
+ -- ({\dx*-3.8400},{\dy*0.3448})
+ -- ({\dx*-3.8200},{\dy*0.3192})
+ -- ({\dx*-3.8000},{\dy*0.2996})
+ -- ({\dx*-3.7800},{\dy*0.2845})
+ -- ({\dx*-3.7600},{\dy*0.2727})
+ -- ({\dx*-3.7400},{\dy*0.2636})
+ -- ({\dx*-3.7200},{\dy*0.2567})
+ -- ({\dx*-3.7000},{\dy*0.2516})
+ -- ({\dx*-3.6800},{\dy*0.2481})
+ -- ({\dx*-3.6600},{\dy*0.2460})
+ -- ({\dx*-3.6400},{\dy*0.2452})
+ -- ({\dx*-3.6200},{\dy*0.2455})
+ -- ({\dx*-3.6000},{\dy*0.2469})
+ -- ({\dx*-3.5800},{\dy*0.2493})
+ -- ({\dx*-3.5600},{\dy*0.2529})
+ -- ({\dx*-3.5400},{\dy*0.2575})
+ -- ({\dx*-3.5200},{\dy*0.2632})
+ -- ({\dx*-3.5000},{\dy*0.2701})
+ -- ({\dx*-3.4800},{\dy*0.2782})
+ -- ({\dx*-3.4600},{\dy*0.2877})
+ -- ({\dx*-3.4400},{\dy*0.2987})
+ -- ({\dx*-3.4200},{\dy*0.3114})
+ -- ({\dx*-3.4000},{\dy*0.3259})
+ -- ({\dx*-3.3800},{\dy*0.3425})
+ -- ({\dx*-3.3600},{\dy*0.3617})
+ -- ({\dx*-3.3400},{\dy*0.3837})
+ -- ({\dx*-3.3200},{\dy*0.4091})
+ -- ({\dx*-3.3000},{\dy*0.4385})
+ -- ({\dx*-3.2800},{\dy*0.4729})
+ -- ({\dx*-3.2600},{\dy*0.5133})
+ -- ({\dx*-3.2400},{\dy*0.5613})
+ -- ({\dx*-3.2200},{\dy*0.6189})
+ -- ({\dx*-3.2000},{\dy*0.6891})
+ -- ({\dx*-3.1800},{\dy*0.7759})
+ -- ({\dx*-3.1600},{\dy*0.8858})
+ -- ({\dx*-3.1400},{\dy*1.0286})
+ -- ({\dx*-3.1200},{\dy*1.2209})
+ -- ({\dx*-3.1000},{\dy*1.4923})
+ -- ({\dx*-3.0800},{\dy*1.9024})
+ -- ({\dx*-3.0600},{\dy*2.5901})
+ -- ({\dx*-3.0400},{\dy*3.9720})
+ -- ({\dx*-3.0200},{\dy*8.1315})
+ -- ({\dx*-3.0050},{\dy*33.1259})}
+\def\gammafive{({\dx*-4.9990},{\dy*-8.3476})
+ -- ({\dx*-4.9800},{\dy*-0.4314})
+ -- ({\dx*-4.9600},{\dy*-0.2236})
+ -- ({\dx*-4.9400},{\dy*-0.1547})
+ -- ({\dx*-4.9200},{\dy*-0.1206})
+ -- ({\dx*-4.9000},{\dy*-0.1004})
+ -- ({\dx*-4.8800},{\dy*-0.0872})
+ -- ({\dx*-4.8600},{\dy*-0.0779})
+ -- ({\dx*-4.8400},{\dy*-0.0712})
+ -- ({\dx*-4.8200},{\dy*-0.0662})
+ -- ({\dx*-4.8000},{\dy*-0.0624})
+ -- ({\dx*-4.7800},{\dy*-0.0595})
+ -- ({\dx*-4.7600},{\dy*-0.0573})
+ -- ({\dx*-4.7400},{\dy*-0.0556})
+ -- ({\dx*-4.7200},{\dy*-0.0544})
+ -- ({\dx*-4.7000},{\dy*-0.0535})
+ -- ({\dx*-4.6800},{\dy*-0.0530})
+ -- ({\dx*-4.6600},{\dy*-0.0528})
+ -- ({\dx*-4.6400},{\dy*-0.0528})
+ -- ({\dx*-4.6200},{\dy*-0.0531})
+ -- ({\dx*-4.6000},{\dy*-0.0537})
+ -- ({\dx*-4.5800},{\dy*-0.0544})
+ -- ({\dx*-4.5600},{\dy*-0.0555})
+ -- ({\dx*-4.5400},{\dy*-0.0567})
+ -- ({\dx*-4.5200},{\dy*-0.0582})
+ -- ({\dx*-4.5000},{\dy*-0.0600})
+ -- ({\dx*-4.4800},{\dy*-0.0621})
+ -- ({\dx*-4.4600},{\dy*-0.0645})
+ -- ({\dx*-4.4400},{\dy*-0.0673})
+ -- ({\dx*-4.4200},{\dy*-0.0704})
+ -- ({\dx*-4.4000},{\dy*-0.0741})
+ -- ({\dx*-4.3800},{\dy*-0.0782})
+ -- ({\dx*-4.3600},{\dy*-0.0830})
+ -- ({\dx*-4.3400},{\dy*-0.0884})
+ -- ({\dx*-4.3200},{\dy*-0.0947})
+ -- ({\dx*-4.3000},{\dy*-0.1020})
+ -- ({\dx*-4.2800},{\dy*-0.1105})
+ -- ({\dx*-4.2600},{\dy*-0.1205})
+ -- ({\dx*-4.2400},{\dy*-0.1324})
+ -- ({\dx*-4.2200},{\dy*-0.1467})
+ -- ({\dx*-4.2000},{\dy*-0.1641})
+ -- ({\dx*-4.1800},{\dy*-0.1856})
+ -- ({\dx*-4.1600},{\dy*-0.2129})
+ -- ({\dx*-4.1400},{\dy*-0.2485})
+ -- ({\dx*-4.1200},{\dy*-0.2963})
+ -- ({\dx*-4.1000},{\dy*-0.3640})
+ -- ({\dx*-4.0800},{\dy*-0.4663})
+ -- ({\dx*-4.0600},{\dy*-0.6380})
+ -- ({\dx*-4.0400},{\dy*-0.9832})
+ -- ({\dx*-4.0200},{\dy*-2.0228})
+ -- ({\dx*-4.0010},{\dy*-41.6040})}
+\def\gammasix{({\dx*-5.9998},{\dy*6.9470})
+ -- ({\dx*-5.9800},{\dy*0.0721})
+ -- ({\dx*-5.9600},{\dy*0.0375})
+ -- ({\dx*-5.9400},{\dy*0.0260})
+ -- ({\dx*-5.9200},{\dy*0.0204})
+ -- ({\dx*-5.9000},{\dy*0.0170})
+ -- ({\dx*-5.8800},{\dy*0.0148})
+ -- ({\dx*-5.8600},{\dy*0.0133})
+ -- ({\dx*-5.8400},{\dy*0.0122})
+ -- ({\dx*-5.8200},{\dy*0.0114})
+ -- ({\dx*-5.8000},{\dy*0.0108})
+ -- ({\dx*-5.7800},{\dy*0.0103})
+ -- ({\dx*-5.7600},{\dy*0.0099})
+ -- ({\dx*-5.7400},{\dy*0.0097})
+ -- ({\dx*-5.7200},{\dy*0.0095})
+ -- ({\dx*-5.7000},{\dy*0.0094})
+ -- ({\dx*-5.6800},{\dy*0.0093})
+ -- ({\dx*-5.6600},{\dy*0.0093})
+ -- ({\dx*-5.6400},{\dy*0.0094})
+ -- ({\dx*-5.6200},{\dy*0.0095})
+ -- ({\dx*-5.6000},{\dy*0.0096})
+ -- ({\dx*-5.5800},{\dy*0.0098})
+ -- ({\dx*-5.5600},{\dy*0.0100})
+ -- ({\dx*-5.5400},{\dy*0.0102})
+ -- ({\dx*-5.5200},{\dy*0.0105})
+ -- ({\dx*-5.5000},{\dy*0.0109})
+ -- ({\dx*-5.4800},{\dy*0.0113})
+ -- ({\dx*-5.4600},{\dy*0.0118})
+ -- ({\dx*-5.4400},{\dy*0.0124})
+ -- ({\dx*-5.4200},{\dy*0.0130})
+ -- ({\dx*-5.4000},{\dy*0.0137})
+ -- ({\dx*-5.3800},{\dy*0.0145})
+ -- ({\dx*-5.3600},{\dy*0.0155})
+ -- ({\dx*-5.3400},{\dy*0.0166})
+ -- ({\dx*-5.3200},{\dy*0.0178})
+ -- ({\dx*-5.3000},{\dy*0.0192})
+ -- ({\dx*-5.2800},{\dy*0.0209})
+ -- ({\dx*-5.2600},{\dy*0.0229})
+ -- ({\dx*-5.2400},{\dy*0.0253})
+ -- ({\dx*-5.2200},{\dy*0.0281})
+ -- ({\dx*-5.2000},{\dy*0.0316})
+ -- ({\dx*-5.1800},{\dy*0.0358})
+ -- ({\dx*-5.1600},{\dy*0.0413})
+ -- ({\dx*-5.1400},{\dy*0.0483})
+ -- ({\dx*-5.1200},{\dy*0.0579})
+ -- ({\dx*-5.1000},{\dy*0.0714})
+ -- ({\dx*-5.0800},{\dy*0.0918})
+ -- ({\dx*-5.0600},{\dy*0.1261})
+ -- ({\dx*-5.0400},{\dy*0.1951})
+ -- ({\dx*-5.0200},{\dy*0.4029})
+ -- ({\dx*-5.0002},{\dy*41.6525})}
+\def\gammasinplus{({\dx*0.0190},{\dy*52.1325})
+ -- ({\dx*0.0200},{\dy*49.5050})
+ -- ({\dx*0.0400},{\dy*24.5863})
+ -- ({\dx*0.0600},{\dy*16.3331})
+ -- ({\dx*0.0800},{\dy*12.2453})
+ -- ({\dx*0.1000},{\dy*9.8225})
+ -- ({\dx*0.1200},{\dy*8.2314})
+ -- ({\dx*0.1400},{\dy*7.1145})
+ -- ({\dx*0.1600},{\dy*6.2930})
+ -- ({\dx*0.1800},{\dy*5.6676})
+ -- ({\dx*0.2000},{\dy*5.1786})
+ -- ({\dx*0.2200},{\dy*4.7879})
+ -- ({\dx*0.2400},{\dy*4.4701})
+ -- ({\dx*0.2600},{\dy*4.2074})
+ -- ({\dx*0.2800},{\dy*3.9874})
+ -- ({\dx*0.3000},{\dy*3.8006})
+ -- ({\dx*0.3200},{\dy*3.6401})
+ -- ({\dx*0.3400},{\dy*3.5005})
+ -- ({\dx*0.3600},{\dy*3.3776})
+ -- ({\dx*0.3800},{\dy*3.2680})
+ -- ({\dx*0.4000},{\dy*3.1692})
+ -- ({\dx*0.4200},{\dy*3.0790})
+ -- ({\dx*0.4400},{\dy*2.9955})
+ -- ({\dx*0.4600},{\dy*2.9173})
+ -- ({\dx*0.4800},{\dy*2.8433})
+ -- ({\dx*0.5000},{\dy*2.7725})
+ -- ({\dx*0.5200},{\dy*2.7039})
+ -- ({\dx*0.5400},{\dy*2.6369})
+ -- ({\dx*0.5600},{\dy*2.5709})
+ -- ({\dx*0.5800},{\dy*2.5055})
+ -- ({\dx*0.6000},{\dy*2.4402})
+ -- ({\dx*0.6200},{\dy*2.3748})
+ -- ({\dx*0.6400},{\dy*2.3090})
+ -- ({\dx*0.6600},{\dy*2.2425})
+ -- ({\dx*0.6800},{\dy*2.1752})
+ -- ({\dx*0.7000},{\dy*2.1071})
+ -- ({\dx*0.7200},{\dy*2.0380})
+ -- ({\dx*0.7400},{\dy*1.9679})
+ -- ({\dx*0.7600},{\dy*1.8969})
+ -- ({\dx*0.7800},{\dy*1.8249})
+ -- ({\dx*0.8000},{\dy*1.7520})
+ -- ({\dx*0.8200},{\dy*1.6783})
+ -- ({\dx*0.8400},{\dy*1.6039})
+ -- ({\dx*0.8600},{\dy*1.5289})
+ -- ({\dx*0.8800},{\dy*1.4534})
+ -- ({\dx*0.9000},{\dy*1.3776})
+ -- ({\dx*0.9200},{\dy*1.3017})
+ -- ({\dx*0.9400},{\dy*1.2258})
+ -- ({\dx*0.9600},{\dy*1.1501})
+ -- ({\dx*0.9800},{\dy*1.0747})
+ -- ({\dx*1.0000},{\dy*1.0000})
+ -- ({\dx*1.0200},{\dy*0.9261})
+ -- ({\dx*1.0400},{\dy*0.8531})
+ -- ({\dx*1.0600},{\dy*0.7814})
+ -- ({\dx*1.0800},{\dy*0.7110})
+ -- ({\dx*1.1000},{\dy*0.6423})
+ -- ({\dx*1.1200},{\dy*0.5755})
+ -- ({\dx*1.1400},{\dy*0.5106})
+ -- ({\dx*1.1600},{\dy*0.4480})
+ -- ({\dx*1.1800},{\dy*0.3879})
+ -- ({\dx*1.2000},{\dy*0.3304})
+ -- ({\dx*1.2200},{\dy*0.2757})
+ -- ({\dx*1.2400},{\dy*0.2240})
+ -- ({\dx*1.2600},{\dy*0.1754})
+ -- ({\dx*1.2800},{\dy*0.1302})
+ -- ({\dx*1.3000},{\dy*0.0885})
+ -- ({\dx*1.3200},{\dy*0.0503})
+ -- ({\dx*1.3400},{\dy*0.0159})
+ -- ({\dx*1.3600},{\dy*-0.0146})
+ -- ({\dx*1.3800},{\dy*-0.0412})
+ -- ({\dx*1.4000},{\dy*-0.0638})
+ -- ({\dx*1.4200},{\dy*-0.0822})
+ -- ({\dx*1.4400},{\dy*-0.0965})
+ -- ({\dx*1.4600},{\dy*-0.1065})
+ -- ({\dx*1.4800},{\dy*-0.1123})
+ -- ({\dx*1.5000},{\dy*-0.1138})
+ -- ({\dx*1.5200},{\dy*-0.1110})
+ -- ({\dx*1.5400},{\dy*-0.1039})
+ -- ({\dx*1.5600},{\dy*-0.0926})
+ -- ({\dx*1.5800},{\dy*-0.0772})
+ -- ({\dx*1.6000},{\dy*-0.0575})
+ -- ({\dx*1.6200},{\dy*-0.0339})
+ -- ({\dx*1.6400},{\dy*-0.0062})
+ -- ({\dx*1.6600},{\dy*0.0254})
+ -- ({\dx*1.6800},{\dy*0.0607})
+ -- ({\dx*1.7000},{\dy*0.0996})
+ -- ({\dx*1.7200},{\dy*0.1421})
+ -- ({\dx*1.7400},{\dy*0.1879})
+ -- ({\dx*1.7600},{\dy*0.2368})
+ -- ({\dx*1.7800},{\dy*0.2888})
+ -- ({\dx*1.8000},{\dy*0.3436})
+ -- ({\dx*1.8200},{\dy*0.4010})
+ -- ({\dx*1.8400},{\dy*0.4609})
+ -- ({\dx*1.8600},{\dy*0.5229})
+ -- ({\dx*1.8800},{\dy*0.5869})
+ -- ({\dx*1.9000},{\dy*0.6527})
+ -- ({\dx*1.9200},{\dy*0.7201})
+ -- ({\dx*1.9400},{\dy*0.7887})
+ -- ({\dx*1.9600},{\dy*0.8584})
+ -- ({\dx*1.9800},{\dy*0.9289})
+ -- ({\dx*2.0000},{\dy*1.0000})
+ -- ({\dx*2.0200},{\dy*1.0714})
+ -- ({\dx*2.0400},{\dy*1.1429})
+ -- ({\dx*2.0600},{\dy*1.2142})
+ -- ({\dx*2.0800},{\dy*1.2852})
+ -- ({\dx*2.1000},{\dy*1.3555})
+ -- ({\dx*2.1200},{\dy*1.4249})
+ -- ({\dx*2.1400},{\dy*1.4933})
+ -- ({\dx*2.1600},{\dy*1.5603})
+ -- ({\dx*2.1800},{\dy*1.6258})
+ -- ({\dx*2.2000},{\dy*1.6896})
+ -- ({\dx*2.2200},{\dy*1.7514})
+ -- ({\dx*2.2400},{\dy*1.8111})
+ -- ({\dx*2.2600},{\dy*1.8685})
+ -- ({\dx*2.2800},{\dy*1.9234})
+ -- ({\dx*2.3000},{\dy*1.9757})
+ -- ({\dx*2.3200},{\dy*2.0253})
+ -- ({\dx*2.3400},{\dy*2.0719})
+ -- ({\dx*2.3600},{\dy*2.1155})
+ -- ({\dx*2.3800},{\dy*2.1560})
+ -- ({\dx*2.4000},{\dy*2.1932})
+ -- ({\dx*2.4200},{\dy*2.2272})
+ -- ({\dx*2.4400},{\dy*2.2578})
+ -- ({\dx*2.4600},{\dy*2.2851})
+ -- ({\dx*2.4800},{\dy*2.3089})
+ -- ({\dx*2.5000},{\dy*2.3293})
+ -- ({\dx*2.5200},{\dy*2.3463})
+ -- ({\dx*2.5400},{\dy*2.3599})
+ -- ({\dx*2.5600},{\dy*2.3701})
+ -- ({\dx*2.5800},{\dy*2.3770})
+ -- ({\dx*2.6000},{\dy*2.3807})
+ -- ({\dx*2.6200},{\dy*2.3812})
+ -- ({\dx*2.6400},{\dy*2.3786})
+ -- ({\dx*2.6600},{\dy*2.3731})
+ -- ({\dx*2.6800},{\dy*2.3647})
+ -- ({\dx*2.7000},{\dy*2.3537})
+ -- ({\dx*2.7200},{\dy*2.3402})
+ -- ({\dx*2.7400},{\dy*2.3242})
+ -- ({\dx*2.7600},{\dy*2.3062})
+ -- ({\dx*2.7800},{\dy*2.2861})
+ -- ({\dx*2.8000},{\dy*2.2643})
+ -- ({\dx*2.8200},{\dy*2.2409})
+ -- ({\dx*2.8400},{\dy*2.2162})
+ -- ({\dx*2.8600},{\dy*2.1903})
+ -- ({\dx*2.8800},{\dy*2.1637})
+ -- ({\dx*2.9000},{\dy*2.1364})
+ -- ({\dx*2.9200},{\dy*2.1087})
+ -- ({\dx*2.9400},{\dy*2.0810})
+ -- ({\dx*2.9600},{\dy*2.0535})
+ -- ({\dx*2.9800},{\dy*2.0264})
+ -- ({\dx*3.0000},{\dy*2.0000})
+ -- ({\dx*3.0200},{\dy*1.9746})
+ -- ({\dx*3.0400},{\dy*1.9505})
+ -- ({\dx*3.0600},{\dy*1.9280})
+ -- ({\dx*3.0800},{\dy*1.9072})
+ -- ({\dx*3.1000},{\dy*1.8886})
+ -- ({\dx*3.1200},{\dy*1.8723})
+ -- ({\dx*3.1400},{\dy*1.8587})
+ -- ({\dx*3.1600},{\dy*1.8480})
+ -- ({\dx*3.1800},{\dy*1.8404})
+ -- ({\dx*3.2000},{\dy*1.8362})
+ -- ({\dx*3.2200},{\dy*1.8356})
+ -- ({\dx*3.2400},{\dy*1.8390})
+ -- ({\dx*3.2600},{\dy*1.8464})
+ -- ({\dx*3.2800},{\dy*1.8581})
+ -- ({\dx*3.3000},{\dy*1.8744})
+ -- ({\dx*3.3200},{\dy*1.8954})
+ -- ({\dx*3.3400},{\dy*1.9213})
+ -- ({\dx*3.3600},{\dy*1.9523})
+ -- ({\dx*3.3800},{\dy*1.9885})
+ -- ({\dx*3.4000},{\dy*2.0301})
+ -- ({\dx*3.4200},{\dy*2.0773})
+ -- ({\dx*3.4400},{\dy*2.1301})
+ -- ({\dx*3.4600},{\dy*2.1886})
+ -- ({\dx*3.4800},{\dy*2.2530})
+ -- ({\dx*3.5000},{\dy*2.3234})
+ -- ({\dx*3.5200},{\dy*2.3997})
+ -- ({\dx*3.5400},{\dy*2.4821})
+ -- ({\dx*3.5600},{\dy*2.5706})
+ -- ({\dx*3.5800},{\dy*2.6652})
+ -- ({\dx*3.6000},{\dy*2.7660})
+ -- ({\dx*3.6200},{\dy*2.8729})
+ -- ({\dx*3.6400},{\dy*2.9859})
+ -- ({\dx*3.6600},{\dy*3.1051})
+ -- ({\dx*3.6800},{\dy*3.2303})
+ -- ({\dx*3.7000},{\dy*3.3616})
+ -- ({\dx*3.7200},{\dy*3.4989})
+ -- ({\dx*3.7400},{\dy*3.6421})
+ -- ({\dx*3.7600},{\dy*3.7911})
+ -- ({\dx*3.7800},{\dy*3.9459})
+ -- ({\dx*3.8000},{\dy*4.1064})
+ -- ({\dx*3.8200},{\dy*4.2724})
+ -- ({\dx*3.8400},{\dy*4.4440})
+ -- ({\dx*3.8600},{\dy*4.6209})
+ -- ({\dx*3.8800},{\dy*4.8030})
+ -- ({\dx*3.9000},{\dy*4.9903})
+ -- ({\dx*3.9200},{\dy*5.1826})
+ -- ({\dx*3.9400},{\dy*5.3799})
+ -- ({\dx*3.9600},{\dy*5.5819})
+ -- ({\dx*3.9800},{\dy*5.7887})
+ -- ({\dx*4.0000},{\dy*6.0000})
+ -- ({\dx*4.0200},{\dy*6.2158})
+ -- ({\dx*4.0400},{\dy*6.4359})
+ -- ({\dx*4.0600},{\dy*6.6603})
+ -- ({\dx*4.0800},{\dy*6.8889})
+ -- ({\dx*4.0810},{\dy*6.9005})}
+\def\gammasinone{({\dx*-0.9810},{\dy*-53.1410})
+ -- ({\dx*-0.9800},{\dy*-50.5140})
+ -- ({\dx*-0.9600},{\dy*-25.6055})
+ -- ({\dx*-0.9400},{\dy*-17.3637})
+ -- ({\dx*-0.9200},{\dy*-13.2884})
+ -- ({\dx*-0.9000},{\dy*-10.8796})
+ -- ({\dx*-0.8800},{\dy*-9.3036})
+ -- ({\dx*-0.8600},{\dy*-8.2033})
+ -- ({\dx*-0.8400},{\dy*-7.3999})
+ -- ({\dx*-0.8200},{\dy*-6.7941})
+ -- ({\dx*-0.8000},{\dy*-6.3263})
+ -- ({\dx*-0.7800},{\dy*-5.9586})
+ -- ({\dx*-0.7600},{\dy*-5.6655})
+ -- ({\dx*-0.7400},{\dy*-5.4296})
+ -- ({\dx*-0.7200},{\dy*-5.2384})
+ -- ({\dx*-0.7000},{\dy*-5.0827})
+ -- ({\dx*-0.6800},{\dy*-4.9557})
+ -- ({\dx*-0.6600},{\dy*-4.8523})
+ -- ({\dx*-0.6400},{\dy*-4.7685})
+ -- ({\dx*-0.6200},{\dy*-4.7012})
+ -- ({\dx*-0.6000},{\dy*-4.6480})
+ -- ({\dx*-0.5800},{\dy*-4.6072})
+ -- ({\dx*-0.5600},{\dy*-4.5773})
+ -- ({\dx*-0.5400},{\dy*-4.5573})
+ -- ({\dx*-0.5200},{\dy*-4.5467})
+ -- ({\dx*-0.5000},{\dy*-4.5449})
+ -- ({\dx*-0.4800},{\dy*-4.5519})
+ -- ({\dx*-0.4600},{\dy*-4.5677})
+ -- ({\dx*-0.4400},{\dy*-4.5928})
+ -- ({\dx*-0.4200},{\dy*-4.6279})
+ -- ({\dx*-0.4000},{\dy*-4.6740})
+ -- ({\dx*-0.3800},{\dy*-4.7325})
+ -- ({\dx*-0.3600},{\dy*-4.8052})
+ -- ({\dx*-0.3400},{\dy*-4.8944})
+ -- ({\dx*-0.3200},{\dy*-5.0033})
+ -- ({\dx*-0.3000},{\dy*-5.1359})
+ -- ({\dx*-0.2800},{\dy*-5.2972})
+ -- ({\dx*-0.2600},{\dy*-5.4942})
+ -- ({\dx*-0.2400},{\dy*-5.7359})
+ -- ({\dx*-0.2200},{\dy*-6.0350})
+ -- ({\dx*-0.2000},{\dy*-6.4089})
+ -- ({\dx*-0.1800},{\dy*-6.8830})
+ -- ({\dx*-0.1600},{\dy*-7.4952})
+ -- ({\dx*-0.1400},{\dy*-8.3052})
+ -- ({\dx*-0.1200},{\dy*-9.4124})
+ -- ({\dx*-0.1000},{\dy*-10.9953})
+ -- ({\dx*-0.0800},{\dy*-13.4114})
+ -- ({\dx*-0.0600},{\dy*-17.4941})
+ -- ({\dx*-0.0400},{\dy*-25.7436})
+ -- ({\dx*-0.0200},{\dy*-50.6602})
+ -- ({\dx*-0.0190},{\dy*-53.2876})}
+\def\gammasintwo{({\dx*-1.9810},{\dy*26.8549})
+ -- ({\dx*-1.9800},{\dy*25.5432})
+ -- ({\dx*-1.9600},{\dy*13.1254})
+ -- ({\dx*-1.9400},{\dy*9.0411})
+ -- ({\dx*-1.9200},{\dy*7.0402})
+ -- ({\dx*-1.9000},{\dy*5.8725})
+ -- ({\dx*-1.8800},{\dy*5.1211})
+ -- ({\dx*-1.8600},{\dy*4.6073})
+ -- ({\dx*-1.8400},{\dy*4.2416})
+ -- ({\dx*-1.8200},{\dy*3.9745})
+ -- ({\dx*-1.8000},{\dy*3.7759})
+ -- ({\dx*-1.7800},{\dy*3.6268})
+ -- ({\dx*-1.7600},{\dy*3.5146})
+ -- ({\dx*-1.7400},{\dy*3.4305})
+ -- ({\dx*-1.7200},{\dy*3.3681})
+ -- ({\dx*-1.7000},{\dy*3.3229})
+ -- ({\dx*-1.6800},{\dy*3.2916})
+ -- ({\dx*-1.6600},{\dy*3.2715})
+ -- ({\dx*-1.6400},{\dy*3.2607})
+ -- ({\dx*-1.6200},{\dy*3.2578})
+ -- ({\dx*-1.6000},{\dy*3.2616})
+ -- ({\dx*-1.5800},{\dy*3.2715})
+ -- ({\dx*-1.5600},{\dy*3.2868})
+ -- ({\dx*-1.5400},{\dy*3.3072})
+ -- ({\dx*-1.5200},{\dy*3.3327})
+ -- ({\dx*-1.5000},{\dy*3.3633})
+ -- ({\dx*-1.4800},{\dy*3.3993})
+ -- ({\dx*-1.4600},{\dy*3.4412})
+ -- ({\dx*-1.4400},{\dy*3.4896})
+ -- ({\dx*-1.4200},{\dy*3.5456})
+ -- ({\dx*-1.4000},{\dy*3.6103})
+ -- ({\dx*-1.3800},{\dy*3.6854})
+ -- ({\dx*-1.3600},{\dy*3.7727})
+ -- ({\dx*-1.3400},{\dy*3.8749})
+ -- ({\dx*-1.3200},{\dy*3.9951})
+ -- ({\dx*-1.3000},{\dy*4.1374})
+ -- ({\dx*-1.2800},{\dy*4.3070})
+ -- ({\dx*-1.2600},{\dy*4.5109})
+ -- ({\dx*-1.2400},{\dy*4.7583})
+ -- ({\dx*-1.2200},{\dy*5.0617})
+ -- ({\dx*-1.2000},{\dy*5.4387})
+ -- ({\dx*-1.1800},{\dy*5.9148})
+ -- ({\dx*-1.1600},{\dy*6.5279})
+ -- ({\dx*-1.1400},{\dy*7.3376})
+ -- ({\dx*-1.1200},{\dy*8.4433})
+ -- ({\dx*-1.1000},{\dy*10.0238})
+ -- ({\dx*-1.0800},{\dy*12.4364})
+ -- ({\dx*-1.0600},{\dy*16.5145})
+ -- ({\dx*-1.0400},{\dy*24.7583})
+ -- ({\dx*-1.0200},{\dy*49.6681})
+ -- ({\dx*-1.0190},{\dy*52.2951})}
+\def\gammasinthree{({\dx*-2.9810},{\dy*-9.0483})
+ -- ({\dx*-2.9800},{\dy*-8.6133})
+ -- ({\dx*-2.9600},{\dy*-4.5173})
+ -- ({\dx*-2.9400},{\dy*-3.1989})
+ -- ({\dx*-2.9200},{\dy*-2.5746})
+ -- ({\dx*-2.9000},{\dy*-2.2274})
+ -- ({\dx*-2.8800},{\dy*-2.0184})
+ -- ({\dx*-2.8600},{\dy*-1.8878})
+ -- ({\dx*-2.8400},{\dy*-1.8057})
+ -- ({\dx*-2.8200},{\dy*-1.7552})
+ -- ({\dx*-2.8000},{\dy*-1.7264})
+ -- ({\dx*-2.7800},{\dy*-1.7127})
+ -- ({\dx*-2.7600},{\dy*-1.7099})
+ -- ({\dx*-2.7400},{\dy*-1.7149})
+ -- ({\dx*-2.7200},{\dy*-1.7255})
+ -- ({\dx*-2.7000},{\dy*-1.7401})
+ -- ({\dx*-2.6800},{\dy*-1.7575})
+ -- ({\dx*-2.6600},{\dy*-1.7768})
+ -- ({\dx*-2.6400},{\dy*-1.7972})
+ -- ({\dx*-2.6200},{\dy*-1.8183})
+ -- ({\dx*-2.6000},{\dy*-1.8397})
+ -- ({\dx*-2.5800},{\dy*-1.8612})
+ -- ({\dx*-2.5600},{\dy*-1.8825})
+ -- ({\dx*-2.5400},{\dy*-1.9036})
+ -- ({\dx*-2.5200},{\dy*-1.9245})
+ -- ({\dx*-2.5000},{\dy*-1.9453})
+ -- ({\dx*-2.4800},{\dy*-1.9663})
+ -- ({\dx*-2.4600},{\dy*-1.9877})
+ -- ({\dx*-2.4400},{\dy*-2.0099})
+ -- ({\dx*-2.4200},{\dy*-2.0335})
+ -- ({\dx*-2.4000},{\dy*-2.0591})
+ -- ({\dx*-2.3800},{\dy*-2.0876})
+ -- ({\dx*-2.3600},{\dy*-2.1200})
+ -- ({\dx*-2.3400},{\dy*-2.1578})
+ -- ({\dx*-2.3200},{\dy*-2.2024})
+ -- ({\dx*-2.3000},{\dy*-2.2561})
+ -- ({\dx*-2.2800},{\dy*-2.3216})
+ -- ({\dx*-2.2600},{\dy*-2.4024})
+ -- ({\dx*-2.2400},{\dy*-2.5032})
+ -- ({\dx*-2.2200},{\dy*-2.6303})
+ -- ({\dx*-2.2000},{\dy*-2.7928})
+ -- ({\dx*-2.1800},{\dy*-3.0032})
+ -- ({\dx*-2.1600},{\dy*-3.2809})
+ -- ({\dx*-2.1400},{\dy*-3.6556})
+ -- ({\dx*-2.1200},{\dy*-4.1772})
+ -- ({\dx*-2.1000},{\dy*-4.9351})
+ -- ({\dx*-2.0800},{\dy*-6.1082})
+ -- ({\dx*-2.0600},{\dy*-8.1132})
+ -- ({\dx*-2.0400},{\dy*-12.2003})
+ -- ({\dx*-2.0200},{\dy*-24.6199})
+ -- ({\dx*-2.0190},{\dy*-25.9316})}
+\def\gammasinfour{({\dx*-3.9950},{\dy*8.4124})
+ -- ({\dx*-3.9800},{\dy*2.2112})
+ -- ({\dx*-3.9600},{\dy*1.2344})
+ -- ({\dx*-3.9400},{\dy*0.9517})
+ -- ({\dx*-3.9200},{\dy*0.8420})
+ -- ({\dx*-3.9000},{\dy*0.8009})
+ -- ({\dx*-3.8800},{\dy*0.7935})
+ -- ({\dx*-3.8600},{\dy*0.8045})
+ -- ({\dx*-3.8400},{\dy*0.8265})
+ -- ({\dx*-3.8200},{\dy*0.8550})
+ -- ({\dx*-3.8000},{\dy*0.8874})
+ -- ({\dx*-3.7800},{\dy*0.9219})
+ -- ({\dx*-3.7600},{\dy*0.9573})
+ -- ({\dx*-3.7400},{\dy*0.9926})
+ -- ({\dx*-3.7200},{\dy*1.0272})
+ -- ({\dx*-3.7000},{\dy*1.0607})
+ -- ({\dx*-3.6800},{\dy*1.0925})
+ -- ({\dx*-3.6600},{\dy*1.1223})
+ -- ({\dx*-3.6400},{\dy*1.1500})
+ -- ({\dx*-3.6200},{\dy*1.1752})
+ -- ({\dx*-3.6000},{\dy*1.1979})
+ -- ({\dx*-3.5800},{\dy*1.2179})
+ -- ({\dx*-3.5600},{\dy*1.2351})
+ -- ({\dx*-3.5400},{\dy*1.2496})
+ -- ({\dx*-3.5200},{\dy*1.2612})
+ -- ({\dx*-3.5000},{\dy*1.2701})
+ -- ({\dx*-3.4800},{\dy*1.2763})
+ -- ({\dx*-3.4600},{\dy*1.2798})
+ -- ({\dx*-3.4400},{\dy*1.2810})
+ -- ({\dx*-3.4200},{\dy*1.2800})
+ -- ({\dx*-3.4000},{\dy*1.2769})
+ -- ({\dx*-3.3800},{\dy*1.2723})
+ -- ({\dx*-3.3600},{\dy*1.2665})
+ -- ({\dx*-3.3400},{\dy*1.2600})
+ -- ({\dx*-3.3200},{\dy*1.2534})
+ -- ({\dx*-3.3000},{\dy*1.2475})
+ -- ({\dx*-3.2800},{\dy*1.2434})
+ -- ({\dx*-3.2600},{\dy*1.2423})
+ -- ({\dx*-3.2400},{\dy*1.2458})
+ -- ({\dx*-3.2200},{\dy*1.2563})
+ -- ({\dx*-3.2000},{\dy*1.2768})
+ -- ({\dx*-3.1800},{\dy*1.3117})
+ -- ({\dx*-3.1600},{\dy*1.3676})
+ -- ({\dx*-3.1400},{\dy*1.4544})
+ -- ({\dx*-3.1200},{\dy*1.5890})
+ -- ({\dx*-3.1000},{\dy*1.8013})
+ -- ({\dx*-3.0800},{\dy*2.1511})
+ -- ({\dx*-3.0600},{\dy*2.7775})
+ -- ({\dx*-3.0400},{\dy*4.0974})
+ -- ({\dx*-3.0200},{\dy*8.1943})
+ -- ({\dx*-3.0050},{\dy*33.1416})}
+\def\gammasinfive{({\dx*-4.9990},{\dy*-8.3507})
+ -- ({\dx*-4.9800},{\dy*-0.4942})
+ -- ({\dx*-4.9600},{\dy*-0.3489})
+ -- ({\dx*-4.9400},{\dy*-0.3421})
+ -- ({\dx*-4.9200},{\dy*-0.3693})
+ -- ({\dx*-4.9000},{\dy*-0.4094})
+ -- ({\dx*-4.8800},{\dy*-0.4553})
+ -- ({\dx*-4.8600},{\dy*-0.5037})
+ -- ({\dx*-4.8400},{\dy*-0.5530})
+ -- ({\dx*-4.8200},{\dy*-0.6021})
+ -- ({\dx*-4.8000},{\dy*-0.6502})
+ -- ({\dx*-4.7800},{\dy*-0.6969})
+ -- ({\dx*-4.7600},{\dy*-0.7418})
+ -- ({\dx*-4.7400},{\dy*-0.7846})
+ -- ({\dx*-4.7200},{\dy*-0.8249})
+ -- ({\dx*-4.7000},{\dy*-0.8626})
+ -- ({\dx*-4.6800},{\dy*-0.8973})
+ -- ({\dx*-4.6600},{\dy*-0.9291})
+ -- ({\dx*-4.6400},{\dy*-0.9577})
+ -- ({\dx*-4.6200},{\dy*-0.9829})
+ -- ({\dx*-4.6000},{\dy*-1.0047})
+ -- ({\dx*-4.5800},{\dy*-1.0230})
+ -- ({\dx*-4.5600},{\dy*-1.0377})
+ -- ({\dx*-4.5400},{\dy*-1.0488})
+ -- ({\dx*-4.5200},{\dy*-1.0563})
+ -- ({\dx*-4.5000},{\dy*-1.0600})
+ -- ({\dx*-4.4800},{\dy*-1.0601})
+ -- ({\dx*-4.4600},{\dy*-1.0566})
+ -- ({\dx*-4.4400},{\dy*-1.0496})
+ -- ({\dx*-4.4200},{\dy*-1.0390})
+ -- ({\dx*-4.4000},{\dy*-1.0251})
+ -- ({\dx*-4.3800},{\dy*-1.0080})
+ -- ({\dx*-4.3600},{\dy*-0.9878})
+ -- ({\dx*-4.3400},{\dy*-0.9647})
+ -- ({\dx*-4.3200},{\dy*-0.9390})
+ -- ({\dx*-4.3000},{\dy*-0.9110})
+ -- ({\dx*-4.2800},{\dy*-0.8810})
+ -- ({\dx*-4.2600},{\dy*-0.8495})
+ -- ({\dx*-4.2400},{\dy*-0.8169})
+ -- ({\dx*-4.2200},{\dy*-0.7841})
+ -- ({\dx*-4.2000},{\dy*-0.7518})
+ -- ({\dx*-4.1800},{\dy*-0.7215})
+ -- ({\dx*-4.1600},{\dy*-0.6947})
+ -- ({\dx*-4.1400},{\dy*-0.6742})
+ -- ({\dx*-4.1200},{\dy*-0.6644})
+ -- ({\dx*-4.1000},{\dy*-0.6730})
+ -- ({\dx*-4.0800},{\dy*-0.7150})
+ -- ({\dx*-4.0600},{\dy*-0.8253})
+ -- ({\dx*-4.0400},{\dy*-1.1085})
+ -- ({\dx*-4.0200},{\dy*-2.0855})
+ -- ({\dx*-4.0010},{\dy*-41.6072})}
+\def\gammasinsix{({\dx*-5.9998},{\dy*6.9477})
+ -- ({\dx*-5.9800},{\dy*0.1349})
+ -- ({\dx*-5.9600},{\dy*0.1629})
+ -- ({\dx*-5.9400},{\dy*0.2134})
+ -- ({\dx*-5.9200},{\dy*0.2691})
+ -- ({\dx*-5.9000},{\dy*0.3260})
+ -- ({\dx*-5.8800},{\dy*0.3829})
+ -- ({\dx*-5.8600},{\dy*0.4391})
+ -- ({\dx*-5.8400},{\dy*0.4940})
+ -- ({\dx*-5.8200},{\dy*0.5472})
+ -- ({\dx*-5.8000},{\dy*0.5985})
+ -- ({\dx*-5.7800},{\dy*0.6477})
+ -- ({\dx*-5.7600},{\dy*0.6945})
+ -- ({\dx*-5.7400},{\dy*0.7387})
+ -- ({\dx*-5.7200},{\dy*0.7800})
+ -- ({\dx*-5.7000},{\dy*0.8184})
+ -- ({\dx*-5.6800},{\dy*0.8537})
+ -- ({\dx*-5.6600},{\dy*0.8856})
+ -- ({\dx*-5.6400},{\dy*0.9142})
+ -- ({\dx*-5.6200},{\dy*0.9392})
+ -- ({\dx*-5.6000},{\dy*0.9606})
+ -- ({\dx*-5.5800},{\dy*0.9783})
+ -- ({\dx*-5.5600},{\dy*0.9923})
+ -- ({\dx*-5.5400},{\dy*1.0024})
+ -- ({\dx*-5.5200},{\dy*1.0086})
+ -- ({\dx*-5.5000},{\dy*1.0109})
+ -- ({\dx*-5.4800},{\dy*1.0094})
+ -- ({\dx*-5.4600},{\dy*1.0039})
+ -- ({\dx*-5.4400},{\dy*0.9947})
+ -- ({\dx*-5.4200},{\dy*0.9816})
+ -- ({\dx*-5.4000},{\dy*0.9648})
+ -- ({\dx*-5.3800},{\dy*0.9443})
+ -- ({\dx*-5.3600},{\dy*0.9203})
+ -- ({\dx*-5.3400},{\dy*0.8929})
+ -- ({\dx*-5.3200},{\dy*0.8621})
+ -- ({\dx*-5.3000},{\dy*0.8283})
+ -- ({\dx*-5.2800},{\dy*0.7914})
+ -- ({\dx*-5.2600},{\dy*0.7519})
+ -- ({\dx*-5.2400},{\dy*0.7098})
+ -- ({\dx*-5.2200},{\dy*0.6655})
+ -- ({\dx*-5.2000},{\dy*0.6193})
+ -- ({\dx*-5.1800},{\dy*0.5717})
+ -- ({\dx*-5.1600},{\dy*0.5230})
+ -- ({\dx*-5.1400},{\dy*0.4741})
+ -- ({\dx*-5.1200},{\dy*0.4260})
+ -- ({\dx*-5.1000},{\dy*0.3804})
+ -- ({\dx*-5.0800},{\dy*0.3405})
+ -- ({\dx*-5.0600},{\dy*0.3135})
+ -- ({\dx*-5.0400},{\dy*0.3204})
+ -- ({\dx*-5.0200},{\dy*0.4657})
+ -- ({\dx*-5.0002},{\dy*41.6531})}
diff --git a/buch/papers/laguerre/images/gammaplot.pdf b/buch/papers/laguerre/images/gammaplot.pdf
new file mode 100644
index 0000000..7c195f2
--- /dev/null
+++ b/buch/papers/laguerre/images/gammaplot.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/gammaplot.tex b/buch/papers/laguerre/images/gammaplot.tex
new file mode 100644
index 0000000..5a68f0a
--- /dev/null
+++ b/buch/papers/laguerre/images/gammaplot.tex
@@ -0,0 +1,73 @@
+%
+% gammaplot.tex -- template for standalon tikz images
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\documentclass[tikz]{standalone}
+\usepackage{amsmath}
+\usepackage{times}
+\usepackage{txfonts}
+\usepackage{pgfplots}
+\usepackage{csvsimple}
+\usetikzlibrary{arrows,intersections,math}
+\input{gammapaths.tex}
+\begin{document}
+\def\skala{1}
+\begin{tikzpicture}[>=latex,thick,scale=\skala]
+
+\definecolor{mainColor}{HTML}{D72864} % OST pink
+
+\draw[->] (-6.1,0) -- (5.3,0) coordinate[label={$z$}];
+\draw[->] (0,-5.1) -- (0,6.4) coordinate[label={right:$\Gamma(z)$}];
+
+\foreach \x in {-1,-2,-3,-4,-5,-6}{
+ \draw (\x,-0.1) -- (\x,0.1);
+ \draw[line width=0.1pt] (\x,-5) -- (\x,6.2);
+}
+\foreach \x in {1,2,3,4,5}{
+ \draw (\x,-0.1) -- (\x,0.1);
+ \node at (\x,0) [below] {$\x$};
+}
+\foreach \y in {-5,-4,-3,-2,-1,1,2,3,4,5,6}{
+ \draw (-0.1,\y) -- (0.1,\y);
+}
+\foreach \y in {1,2,3,4,5,6}{
+ \node at (0,\y) [left] {$\y$};
+}
+\foreach \y in {-1,-2,-3,-4,-5}{
+ \node at (0,\y) [right] {$\y$};
+}
+\foreach \x in {-1,-3,-5}{
+ \node at (\x,0) [below left] {$\x$};
+}
+\foreach \x in {-2,-4,-6}{
+ \node at (\x,0) [above left] {$\x$};
+}
+
+\def\dx{1}
+\def\dy{1}
+
+\begin{scope}
+\clip (-6.1,-5) rectangle (4.3,6.2);
+
+% \draw[color=darkgreen,line width=1.4pt] \gammasinplus;
+% \draw[color=darkgreen,line width=1.4pt] \gammasinone;
+% \draw[color=darkgreen,line width=1.4pt] \gammasintwo;
+% \draw[color=darkgreen,line width=1.4pt] \gammasinthree;
+% \draw[color=darkgreen,line width=1.4pt] \gammasinfour;
+% \draw[color=darkgreen,line width=1.4pt] \gammasinfive;
+% \draw[color=darkgreen,line width=1.4pt] \gammasinsix;
+
+\draw[color=mainColor,line width=1.4pt] \gammaplus;
+\draw[color=mainColor,line width=1.4pt] \gammaone;
+\draw[color=mainColor,line width=1.4pt] \gammatwo;
+\draw[color=mainColor,line width=1.4pt] \gammathree;
+\draw[color=mainColor,line width=1.4pt] \gammafour;
+\draw[color=mainColor,line width=1.4pt] \gammafive;
+\draw[color=mainColor,line width=1.4pt] \gammasix;
+
+\end{scope}
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/buch/papers/laguerre/images/integrand.pdf b/buch/papers/laguerre/images/integrand.pdf
new file mode 100644
index 0000000..76be412
--- /dev/null
+++ b/buch/papers/laguerre/images/integrand.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/integrand_exp.pdf b/buch/papers/laguerre/images/integrand_exp.pdf
new file mode 100644
index 0000000..5fe7a7b
--- /dev/null
+++ b/buch/papers/laguerre/images/integrand_exp.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/laguerre_poly.pdf b/buch/papers/laguerre/images/laguerre_poly.pdf
new file mode 100644
index 0000000..21278f5
--- /dev/null
+++ b/buch/papers/laguerre/images/laguerre_poly.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/laguerre_polynomes.pdf b/buch/papers/laguerre/images/laguerre_polynomes.pdf
deleted file mode 100644
index 3976bc7..0000000
--- a/buch/papers/laguerre/images/laguerre_polynomes.pdf
+++ /dev/null
Binary files differ
diff --git a/buch/papers/laguerre/images/rel_error_complex.pdf b/buch/papers/laguerre/images/rel_error_complex.pdf
new file mode 100644
index 0000000..c7bb37a
--- /dev/null
+++ b/buch/papers/laguerre/images/rel_error_complex.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/rel_error_mirror.pdf b/buch/papers/laguerre/images/rel_error_mirror.pdf
new file mode 100644
index 0000000..8a27d41
--- /dev/null
+++ b/buch/papers/laguerre/images/rel_error_mirror.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/rel_error_range.pdf b/buch/papers/laguerre/images/rel_error_range.pdf
new file mode 100644
index 0000000..bb8a2d7
--- /dev/null
+++ b/buch/papers/laguerre/images/rel_error_range.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/rel_error_shifted.pdf b/buch/papers/laguerre/images/rel_error_shifted.pdf
new file mode 100644
index 0000000..b7e72dc
--- /dev/null
+++ b/buch/papers/laguerre/images/rel_error_shifted.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/rel_error_simple.pdf b/buch/papers/laguerre/images/rel_error_simple.pdf
new file mode 100644
index 0000000..3212e42
--- /dev/null
+++ b/buch/papers/laguerre/images/rel_error_simple.pdf
Binary files differ
diff --git a/buch/papers/laguerre/images/targets.pdf b/buch/papers/laguerre/images/targets.pdf
new file mode 100644
index 0000000..9514a6d
--- /dev/null
+++ b/buch/papers/laguerre/images/targets.pdf
Binary files differ
diff --git a/buch/papers/laguerre/main.tex b/buch/papers/laguerre/main.tex
index 00e3b43..57a6560 100644
--- a/buch/papers/laguerre/main.tex
+++ b/buch/papers/laguerre/main.tex
@@ -8,7 +8,25 @@
\begin{refsection}
\chapterauthor{Patrik Müller}
-{\large \color{red} TODO: Einleitung}
+{\parindent0pt Die} Laguerre\--Polynome,
+benannt nach Edmond Laguerre (1834 - 1886),
+sind Lösungen der ebenfalls nach Laguerre benannten Differentialgleichung.
+Laguerre entdeckte diese Polynome, als er Approximations\-methoden
+für das Integral
+% $\int_0^\infty \exp(-x) / x \, dx $
+\begin{align*}
+\int_0^\infty \frac{e^{-x}}{x} \, dx
+\end{align*}
+suchte.
+Darum möchten wir uns in diesem Kapitel,
+ganz im Sinne des Entdeckers,
+den Laguerre-Polynomen für Approximationen von Integralen mit
+exponentiell-abfallenden Funktionen widmen.
+Namentlich werden wir versuchen, mittels Laguerre-Polynomen und
+der Gauss-Quadratur eine geeignete Approximation für die Gamma-Funktion zu finden.
+
+Laguerre-Polynome tauchen zudem auch in der Quantenmechanik im radialen Anteil
+der Lösung für die Schrödinger-Gleichung eines Wasserstoffatoms auf.
\input{papers/laguerre/definition}
\input{papers/laguerre/eigenschaften}
diff --git a/buch/papers/laguerre/packages.tex b/buch/papers/laguerre/packages.tex
index 4ebc172..a80d091 100644
--- a/buch/papers/laguerre/packages.tex
+++ b/buch/papers/laguerre/packages.tex
@@ -6,4 +6,4 @@
% if your paper needs special packages, add package commands as in the
% following example
-\usepackage{derivative}
+\DeclareMathOperator{\real}{Re} \ No newline at end of file
diff --git a/buch/papers/laguerre/presentation/presentation.pdf b/buch/papers/laguerre/presentation/presentation.pdf
new file mode 100644
index 0000000..3d00de3
--- /dev/null
+++ b/buch/papers/laguerre/presentation/presentation.pdf
Binary files differ
diff --git a/buch/papers/laguerre/presentation/presentation.tex b/buch/papers/laguerre/presentation/presentation.tex
new file mode 100644
index 0000000..3db69f5
--- /dev/null
+++ b/buch/papers/laguerre/presentation/presentation.tex
@@ -0,0 +1,134 @@
+\documentclass[ngerman, aspectratio=169, xcolor={rgb}]{beamer}
+
+% style
+\mode<presentation>{
+ \usetheme{Frankfurt}
+}
+%packages
+\usepackage[utf8]{inputenc}\DeclareUnicodeCharacter{2212}{-}
+\usepackage[ngerman]{babel}
+\usepackage{graphicx}
+\usepackage{array}
+
+\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
+\usepackage{ragged2e}
+
+\usepackage{bm} % bold math
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{mathtools}
+\usepackage{amsmath}
+\usepackage{multirow} % multi row in tables
+\usepackage{booktabs} %toprule midrule bottomrue in tables
+\usepackage{scrextend}
+\usepackage{textgreek}
+\usepackage[rgb]{xcolor}
+
+\usepackage{ marvosym } % \Lightning
+
+\usepackage{multimedia} % embedded videos
+
+\usepackage{tikz}
+\usepackage{pgf}
+\usepackage{pgfplots}
+
+\usepackage{algorithmic}
+
+%citations
+\usepackage[style=verbose,backend=biber]{biblatex}
+\addbibresource{references.bib}
+
+
+%math font
+\usefonttheme[onlymath]{serif}
+
+%Beamer Template modifications
+%\definecolor{mainColor}{HTML}{0065A3} % HSR blue
+\definecolor{mainColor}{HTML}{D72864} % OST pink
+\definecolor{invColor}{HTML}{28d79b} % OST pink
+\definecolor{dgreen}{HTML}{38ad36} % Dark green
+
+%\definecolor{mainColor}{HTML}{000000} % HSR blue
+\setbeamercolor{palette primary}{bg=white,fg=mainColor}
+\setbeamercolor{palette secondary}{bg=orange,fg=mainColor}
+\setbeamercolor{palette tertiary}{bg=yellow,fg=red}
+\setbeamercolor{palette quaternary}{bg=mainColor,fg=white} %bg = Top bar, fg = active top bar topic
+\setbeamercolor{structure}{fg=black} % itemize, enumerate, etc (bullet points)
+\setbeamercolor{section in toc}{fg=black} % TOC sections
+\setbeamertemplate{section in toc}[sections numbered]
+\setbeamertemplate{subsection in toc}{%
+ \hspace{1.2em}{$\bullet$}~\inserttocsubsection\par}
+
+\setbeamertemplate{itemize items}[circle]
+\setbeamertemplate{description item}[circle]
+\setbeamertemplate{title page}[default][colsep=-4bp,rounded=true]
+\beamertemplatenavigationsymbolsempty
+
+\setbeamercolor{footline}{fg=gray}
+\setbeamertemplate{footline}{%
+ \hfill\usebeamertemplate***{navigation symbols}
+ \hspace{0.5cm}
+ \insertframenumber{}\hspace{0.2cm}\vspace{0.2cm}
+}
+
+\usepackage{caption}
+\captionsetup{labelformat=empty}
+
+%Title Page
+\title{Laguerre-Polynome}
+\subtitle{Anwendung: Approximation der Gamma-Funktion}
+\author{Patrik Müller}
+% \institute{OST Ostschweizer Fachhochschule}
+% \institute{\includegraphics[scale=0.3]{../img/ost_logo.png}}
+\date{\today}
+
+\input{../packages.tex}
+
+\newcommand*{\QED}{\hfill\ensuremath{\blacksquare}}%
+
+\newcommand*{\HL}{\textcolor{mainColor}}
+\newcommand*{\RD}{\textcolor{red}}
+\newcommand*{\BL}{\textcolor{blue}}
+\newcommand*{\GN}{\textcolor{dgreen}}
+
+\definecolor{darkgreen}{rgb}{0,0.6,0}
+
+
+\makeatletter
+\newcount\my@repeat@count
+\newcommand{\myrepeat}[2]{%
+ \begingroup
+ \my@repeat@count=\z@
+ \@whilenum\my@repeat@count<#1\do{#2\advance\my@repeat@count\@ne}%
+ \endgroup
+}
+\makeatother
+
+\usetikzlibrary{automata,arrows,positioning,calc,shapes.geometric, fadings}
+
+\begin{document}
+
+\begin{frame}
+ \titlepage
+\end{frame}
+
+\begin{frame}{Inhaltsverzeichnis}
+ \tableofcontents
+\end{frame}
+
+\input{sections/laguerre}
+
+\input{sections/gaussquad}
+
+\input{sections/gamma}
+
+\input{sections/gamma_approx}
+
+\appendix
+\begin{frame}
+ % \centering
+ % \Large
+ % \textbf{Vielen Dank für die Aufmerksamkeit}
+\end{frame}
+
+\end{document}
diff --git a/buch/papers/laguerre/presentation/sections/gamma.tex b/buch/papers/laguerre/presentation/sections/gamma.tex
new file mode 100644
index 0000000..7dca39b
--- /dev/null
+++ b/buch/papers/laguerre/presentation/sections/gamma.tex
@@ -0,0 +1,51 @@
+\section{Gamma-Funktion}
+
+\begin{frame}{Gamma-Funktion}
+\begin{columns}
+
+\begin{column}{0.55\textwidth}
+\begin{figure}[h]
+\vspace{-16pt}
+\centering
+% \scalebox{0.51}{\input{../images/gammaplot.pdf}}
+\includegraphics[width=1\textwidth]{../images/gammaplot.pdf}
+% \caption{Gamma-Funktion}
+\end{figure}
+\end{column}
+
+\begin{column}{0.45\textwidth}
+Verallgemeinerung der Fakultät
+\begin{align*}
+\Gamma(n) = (n-1)!
+\end{align*}
+
+Integralformel
+\begin{align*}
+\Gamma(z)
+=
+\int_0^\infty x^{z-1} e^{-x} \, dx
+,\quad
+\operatorname{Re} z > 0
+\end{align*}
+
+Funktionalgleichung
+\begin{align*}
+z \Gamma(z)
+=
+\Gamma(z + 1)
+\end{align*}
+
+Reflektionsformel
+\begin{align*}
+\Gamma(z) \Gamma(1 - z)
+=
+\frac{\pi}{\sin \pi z}
+, \quad
+\text{für }
+z \notin \mathbb{Z}
+\end{align*}
+
+\end{column}
+\end{columns}
+
+\end{frame} \ No newline at end of file
diff --git a/buch/papers/laguerre/presentation/sections/gamma_approx.tex b/buch/papers/laguerre/presentation/sections/gamma_approx.tex
new file mode 100644
index 0000000..ecd02ab
--- /dev/null
+++ b/buch/papers/laguerre/presentation/sections/gamma_approx.tex
@@ -0,0 +1,201 @@
+\section{Approximieren der Gamma-Funktion}
+
+\begin{frame}{Anwenden der Gauss-Laguerre-Quadratur auf $\Gamma(z)$}
+
+\begin{align*}
+\Gamma(z)
+ & =
+\int_0^\infty x^{z-1} e^{-x} \, dx
+\uncover<2->{
+\approx
+\sum_{i=1}^{n} f(x_i) A_i
+}
+\uncover<3->{
+=
+\sum_{i=1}^{n} x^{z-1} A_i
+}
+\\\\
+\uncover<4->{
+ & \text{wobei }
+A_i = \frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2}
+\text{ und $x_i$ die Nullstellen von $L_n(x)$}
+}
+\end{align*}
+
+\end{frame}
+
+\begin{frame}{Fehlerabschätzung}
+\begin{align*}
+R_n(\xi)
+ & =
+\frac{(n!)^2}{(2n)!} f^{(2n)}(\xi)
+\\
+ & =
+(z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z - 2n - 1}
+,\quad
+0 < \xi < \infty
+\end{align*}
+
+% \textbf{Probleme:}
+\begin{itemize}
+\item Funktion ist unbeschränkt
+\item Maximum von $R_n$ gibt oberes Limit des Fehlers an
+\uncover<2->{\item[$\Rightarrow$] Schwierig ein Maximum von $R_n(\xi)$ zu finden}
+\end{itemize}
+\end{frame}
+
+\begin{frame}{Einfacher Ansatz}
+
+\begin{figure}[h]
+\centering
+% \scalebox{0.91}{\input{../images/rel_error_simple.pgf}}
+% \resizebox{!}{0.72\textheight}{\input{../images/rel_error_simple.pgf}}
+\includegraphics[width=0.77\textwidth]{../images/rel_error_simple.pdf}
+\caption{Relativer Fehler des einfachen Ansatzes für verschiedene reele Werte
+von $z$ und Grade $n$ der Laguerre-Polynome}
+\end{figure}
+
+\end{frame}
+
+\begin{frame}{Wieso sind die Resultate so schlecht?}
+
+\textbf{Beobachtungen}
+\begin{itemize}
+\item Wenn $z \in \mathbb{Z}$ relativer Fehler $\rightarrow 0$
+\item Gewisse Periodizität zu erkennen
+\item Für grosse und kleine $z$ ergibt sich ein schlechter relativer Fehler
+\item Es gibt Intervalle $[a,a+1]$ mit minimalem relativem Fehler
+\item $a$ ist abhängig von $n$
+\end{itemize}
+
+\uncover<2->{
+\textbf{Ursache?}
+\begin{itemize}
+\item Vermutung: Integrand ist problematisch
+}
+\uncover<3->{
+\item[$\Rightarrow$] Analysieren von $f(x)$ und dem Integranden
+}
+\end{itemize}
+\end{frame}
+
+\begin{frame}{$f(x) = x^z$}
+\begin{figure}[h]
+\centering
+% \scalebox{0.91}{\input{../images/integrand.pgf}}
+\includegraphics[width=0.8\textwidth]{../images/integrand.pdf}
+% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$}
+\end{figure}
+\end{frame}
+
+\begin{frame}{Integrand $x^z e^{-x}$}
+\begin{figure}[h]
+\centering
+% \scalebox{0.91}{\input{../images/integrand_exp.pgf}}
+\includegraphics[width=0.8\textwidth]{../images/integrand_exp.pdf}
+% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$}
+\end{figure}
+\end{frame}
+
+\begin{frame}{Neuer Ansatz?}
+
+\textbf{Vermutung}
+\begin{itemize}
+\item Es gibt Intervalle $[a(n), a(n)+1]$ in denen der relative Fehler minimal
+ist
+\item $a(n) > 0$
+\end{itemize}
+
+\uncover<2->{
+\textbf{Idee}
+\begin{itemize}
+\item[$\Rightarrow$] Berechnen von $\Gamma(z)$ im geeigneten Intervall und dann
+mit Funktionalgleichung zurückverschieben
+\end{itemize}
+}
+
+\uncover<3->{
+\textbf{Wie finden wir $\boldsymbol{a(n)}$?}
+\begin{itemize}
+\item Minimieren des Fehlerterms mit zusätzlichem Verschiebungsterm
+}
+\uncover<4->{$\Rightarrow$ Schwierig das Maximum des Fehlerterms zu bestimmen}
+\uncover<5->{\item Empirisch $a(n)$ bestimmen}
+\uncover<6->{$\Rightarrow$ Sinnvoll,
+da Gauss-Quadratur nur für kleine $n$ praktischen Nutzen hat}
+\end{itemize}
+\end{frame}
+
+\begin{frame}{Verschiebungsterm}
+\begin{columns}
+\begin{column}{0.625\textwidth}
+\begin{figure}[h]
+\centering
+\includegraphics[width=1\textwidth]{../images/targets.pdf}
+\caption{Optimaler Verschiebungsterm $m^*$ in Abhängigkeit von $z$ und $n$}
+\end{figure}
+\end{column}
+\begin{column}{0.375\textwidth}
+\begin{align*}
+\Gamma(z)
+\approx
+\frac{1}{(z-m)_{m}} \sum_{i=1}^{n} x_i^{z + m - 1} A_i
+\end{align*}
+\end{column}
+\end{columns}
+\end{frame}
+
+\begin{frame}{Schätzen von $m^*$}
+\begin{columns}
+\begin{column}{0.65\textwidth}
+\begin{figure}
+\centering
+\vspace{-12pt}
+% \scalebox{0.7}{\input{../images/estimates.pgf}}
+\includegraphics[width=1.0\textwidth]{../images/estimates.pdf}
+% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$}
+\end{figure}
+\end{column}
+\begin{column}{0.34\textwidth}
+\begin{align*}
+\hat{m}
+&=
+\alpha n + \beta
+\\
+&\approx
+1.34093 n + 0.854093
+\\
+m^*
+&=
+\lceil \hat{m} - \operatorname{Re}z \rceil
+\end{align*}
+\end{column}
+\end{columns}
+
+\end{frame}
+
+\begin{frame}{}
+\begin{figure}[h]
+\centering
+% \scalebox{0.6}{\input{../images/rel_error_shifted.pgf}}
+\includegraphics{../images/rel_error_shifted.pdf}
+\caption{Relativer Fehler mit $n=8$, unterschiedlichen Verschiebungstermen $m$ und $z\in(0, 1)$}
+\end{figure}
+\end{frame}
+
+\begin{frame}{}
+\begin{figure}[h]
+\centering
+% \scalebox{0.6}{\input{../images/rel_error_range.pgf}}
+\includegraphics{../images/rel_error_range.pdf}
+\caption{Relativer Fehler mit $n=8$, Verschiebungsterm $m^*$ und $z\in(-5, 5)$}
+\end{figure}
+\end{frame}
+
+\begin{frame}{Vergleich mit Lanczos-Methode}
+Maximaler relativer Fehler für $n=6$
+\begin{itemize}
+ \item Lanczos-Methode $< 10^{-12}$
+ \item Unsere Methode $\approx 10^{-6}$
+\end{itemize}
+\end{frame} \ No newline at end of file
diff --git a/buch/papers/laguerre/presentation/sections/gaussquad.tex b/buch/papers/laguerre/presentation/sections/gaussquad.tex
new file mode 100644
index 0000000..4d973b8
--- /dev/null
+++ b/buch/papers/laguerre/presentation/sections/gaussquad.tex
@@ -0,0 +1,67 @@
+\section{Gauss-Quadratur}
+
+\begin{frame}{Gauss-Quadratur}
+\textbf{Idee}
+\begin{itemize}[<+->]
+\item Polynome können viele Funktionen approximieren
+\item Wenn Verfahren gut für Polynome funktioniert,
+sollte es auch für andere Funktionen funktionieren
+\item Integrieren eines Interpolationspolynom
+\item Interpolationspolynom ist durch Funktionswerte $f(x_i)$ bestimmt
+$\Rightarrow$ Integral kann durch Funktionswerte berechnet werden
+\item Evaluation der Funktionswerte an geeigneten Stellen
+\end{itemize}
+\end{frame}
+
+\begin{frame}{Gauss-Quadratur}
+\begin{align*}
+\int_{-1}^{1} f(x) \, dx
+\approx
+\sum_{i=1}^n f(x_i) A_i
+\end{align*}
+
+\begin{itemize}[<+->]
+\item Exakt für Polynome mit Grad $2n-1$
+\item Interpolationspolynome müssen orthogonal sein
+\item Stützstellen $x_i$ sind Nullstellen des Polynoms
+\item Fehler:
+\begin{align*}
+E
+=
+\frac{f^{(2n)}(\xi)}{(2n)!} \int_{-1}^{1} l(x)^2 \, dx
+,\quad
+\text{wobei }
+l(x) = \prod_{i=1}^n (x-x_i)
+\end{align*}
+\end{itemize}
+\end{frame}
+
+\begin{frame}{Gauss-Laguerre-Quadratur}
+\begin{itemize}[<+->]
+\item Erweiterung des Integrationsintervall von $[-1, 1]$ auf $(a, b)$
+\item Hinzufügen einer Gewichtsfunktion
+\item Bei uneigentlichen Integralen muss Gewichtsfunktion schneller als jedes
+Integrationspolynom gegen $0$ gehen
+\item[$\Rightarrow$] Für Laguerre-Polynome haben wir den Definitionsbereich
+$(0, \infty)$ und die Gewichtsfunktion $w(x) = e^{-x}$
+\begin{align*}
+\int_0^\infty & f(x) e^{-x} \, dx
+\approx
+\sum_{i=1}^n f(x_i) A_i
+\\
+ & \text{wobei }
+A_i = \frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2}
+\text{ und $x_i$ die Nullstellen von $L_n(x)$}
+\end{align*}
+\end{itemize}
+\end{frame}
+
+\begin{frame}{Fehler der Gauss-Laguerre-Quadratur}
+\begin{align*}
+R_n
+=
+\frac{(n!)^2}{(2n)!} f^{(2n)}(\xi)
+,\quad
+0 < \xi < \infty
+\end{align*}
+\end{frame} \ No newline at end of file
diff --git a/buch/papers/laguerre/presentation/sections/laguerre.tex b/buch/papers/laguerre/presentation/sections/laguerre.tex
new file mode 100644
index 0000000..f99214e
--- /dev/null
+++ b/buch/papers/laguerre/presentation/sections/laguerre.tex
@@ -0,0 +1,91 @@
+\section{Laguerre-Polynome}
+
+\begin{frame}{Laguerre-Differentialgleichung}
+
+\begin{itemize}
+\item Benannt nach Edmond Nicolas Laguerre (1834-1886)
+\item Aus Artikel von 1879,
+in dem er $\int_0^\infty \exp(-x)/x \, dx$ analysierte
+\end{itemize}
+
+\begin{align*}
+x y''(x) + (1 - x) y'(x) + n y(x)
+ & =
+0
+, \quad
+n \in \mathbb{N}_0
+, \quad
+x \in \mathbb{R}
+\end{align*}
+
+\end{frame}
+
+\begin{frame}{Lösen der Differentialgleichung}
+
+\begin{align*}
+x y''(x) + (1 - x) y'(x) + n y(x)
+ & =
+0
+\\
+\end{align*}
+
+\uncover<2->{
+\centering
+\begin{tikzpicture}[remember picture,overlay]
+%% use here too
+\path[draw=mainColor, very thick,->](0, 1.1) to
+node[anchor=west]{Potenzreihenansatz} (0, -0.8);
+\end{tikzpicture}
+}
+
+\begin{align*}
+\uncover<3->{
+L_n(x)
+ & =
+\sum_{k=0}^{n} \frac{(-1)^k}{k!} \binom{n}{k} x^k
+}
+\end{align*}
+\uncover<4->{
+\begin{itemize}
+ \item Die Lösungen der DGL sind die Laguerre-Polynome
+\end{itemize}
+}
+\end{frame}
+
+\begin{frame}
+\begin{figure}[h]
+\centering
+% \resizebox{0.74\textwidth}{!}{\input{../images/laguerre_poly.pgf}}
+\includegraphics[width=0.7\textwidth]{../images/laguerre_poly.pdf}
+\caption{Laguerre-Polynome vom Grad $0$ bis $7$}
+\end{figure}
+\end{frame}
+
+\begin{frame}{Orthogonalität}
+\begin{itemize}[<+->]
+\item Beweis: Umformen in Sturm-Liouville-Problem (siehe Paper)
+\begin{alignat*}{5}
+((p(x) &y'(x)))' + q(x) &y(x)
+&=
+\lambda &w(x) &y(x)
+\\
+((x e^{-x} &y'(x)))' + 0 &y(x)
+&=
+n &e^{-x} &y(x)
+\end{alignat*}
+\item Definitionsbereich $(0, \infty)$
+\item Gewichtsfunktion $w(x) = e^{-x}$
+\end{itemize}
+
+\uncover<4->{
+\begin{align*}
+\int_0^\infty e^{-x} L_n(x) L_m(x) \, dx
+=
+0
+,\quad
+n \neq m
+,\quad
+n, m \in \mathbb{N}
+\end{align*}
+}
+\end{frame} \ No newline at end of file
diff --git a/buch/papers/laguerre/quadratur.tex b/buch/papers/laguerre/quadratur.tex
index 60fad7f..a494362 100644
--- a/buch/papers/laguerre/quadratur.tex
+++ b/buch/papers/laguerre/quadratur.tex
@@ -5,35 +5,70 @@
%
\section{Gauss-Quadratur
\label{laguerre:section:quadratur}}
- {\large \color{red} TODO: Einleitung und kurze Beschreibung Gauss-Quadratur}
+Die Gauss-Quadratur ist ein numerisches Integrationsverfahren,
+welches die Eigenschaften von orthogonalen Polynomen verwendet.
+Herleitungen und Analysen der Gauss-Quadratur können im
+Abschnitt~\ref{buch:orthogonal: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 funktioniert,
+sollte es auch für andere Funktionen angemessene Resultate liefern.
+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 einer gewichteten Summe der Form
\begin{align}
-\int_a^b f(x) w(x)
+\int_a^b f(x) w(x) \, dx
\approx
-\sum_{i=1}^N f(x_i) A_i
+\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 wiefolgt 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,
+kann dies behoben werden.
+Für unseren Fall gilt $a = 0$.
+Das Integral eines Polynomes in diesem Intervall ist immer divergent.
+Darum müssen wir das Polynom mit einer Funktion multiplizieren,
+die schneller als jedes Polynom gegen $0$ geht,
+damit das Integral immer noch konvergiert.
+Die Laguerre-Polynome $L_n$ schaffen hier Abhilfe,
+da ihre Gewichtsfunktion $w(x) = 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}$.
+Die Gleichung~\eqref{laguerre:gaussquadratur} lässt sich wie folgt
+umformulieren:
\begin{align}
\int_{0}^{\infty} f(x) e^{-x} dx
\approx
-\sum_{i=1}^{N} f(x_i) A_i
+\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.
+Für das Laguerre-Polynom $L_n$ müssen demnach dessen Nullstellen $x_i$ und
+als Gewichte $A_i$ die Integrale $l_i(x)e^{-x}$ verwendet werden.
Dabei sind
\begin{align*}
l_i(x_j)
@@ -41,39 +76,112 @@ l_i(x_j)
\delta_{ij}
=
\begin{cases}
-1 & i=j \\
-0 & \text{sonst.}
+1 & i=j \\
+0 & \text{sonst}
\end{cases}
+% .
\end{align*}
-Laut \cite{abramowitz+stegun} sind die Gewichte also
-\begin{align}
+die Lagrangschen Interpolationspolynome.
+Laut \cite{laguerre:hildebrand2013introduction} können die Gewichte 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 Laguerre-Koeffizienten 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,
+vereinfacht sich der Term zu
+\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 sich
+\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
+\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
+und \cite{laguerre: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{lagurre:lag_error}
+\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 6956ade..d21009b 100644
--- a/buch/papers/laguerre/references.bib
+++ b/buch/papers/laguerre/references.bib
@@ -3,8 +3,17 @@
%
% (c) 2020 Autor, Hochschule Rapperswil
%
+@book{laguerre:hildebrand2013introduction,
+ title={Introduction to Numerical Analysis: Second Edition},
+ author={Hildebrand, F.B.},
+ isbn={9780486318554},
+ series={Dover Books on Mathematics},
+ year={2013},
+ publisher={Dover Publications},
+ pages = {389}
+}
-@book{abramowitz+stegun,
+@book{laguerre:abramowitz+stegun,
added-at = {2008-06-25T06:25:58.000+0200},
address = {New York},
author = {Abramowitz, Milton and Stegun, Irene A.},
@@ -19,4 +28,23 @@
timestamp = {2008-06-25T06:25:58.000+0200},
title = {Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables},
year = 1972
+}
+
+@article{laguerre:Cassity1965AbcissasCA,
+ title={Abcissas, coefficients, and error term for the generalized Gauss-Laguerre quadrature formula using the zero ordinate},
+ author={C. Ronald Cassity},
+ journal={Mathematics of Computation},
+ year={1965},
+ volume={19},
+ pages={287-296}
+}
+
+@online{laguerre:lanczos,
+ title = {Lanczos Approximation},
+ author={Eric W. Weisstein},
+ url = {https://mathworld.wolfram.com/LanczosApproximation.html},
+ date = {2022-07-18},
+ year = {2022},
+ month = {7},
+ day = {18}
} \ No newline at end of file
diff --git a/buch/papers/laguerre/scripts/estimates.py b/buch/papers/laguerre/scripts/estimates.py
new file mode 100644
index 0000000..21551f3
--- /dev/null
+++ b/buch/papers/laguerre/scripts/estimates.py
@@ -0,0 +1,49 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ import gamma_approx as ga
+ import targets
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ N = 200
+ ns = np.arange(2, 13)
+ step = 1 / (N - 1)
+ x = np.linspace(step, 1 - step, N + 1)
+
+ bests = targets.find_best_loc(N, ns=ns)
+ mean_m = np.mean(bests, -1)
+
+ intercept, bias = np.polyfit(ns, mean_m, 1)
+ fig, axs = plt.subplots(
+ 2, num=1, sharex=True, clear=True, constrained_layout=True, figsize=(4.5, 3.6)
+ )
+ xl = np.array([ns[0] - 0.5, ns[-1] + 0.5])
+ axs[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$")
+ axs[0].plot(ns, mean_m, "x", label=r"$\overline{m}$")
+ axs[1].plot(
+ ns, ((intercept * ns + bias) - mean_m), "-x", label=r"$\hat{m} - \overline{m}$"
+ )
+ axs[0].set_xlim(*xl)
+ axs[0].set_xticks(ns)
+ axs[0].set_yticks(np.arange(np.floor(mean_m[0]), np.ceil(mean_m[-1]) + 0.1, 2))
+ # axs[0].set_title("Schätzung von Mittelwert")
+ # axs[1].set_title("Fehler")
+ axs[-1].set_xlabel(r"$n$")
+ for ax in axs:
+ ax.grid(1)
+ ax.legend()
+ # fig.savefig(f"{ga.img_path}/estimates.pgf")
+ fig.savefig(f"{ga.img_path}/estimates.pdf")
+
+ print(f"Intercept={intercept:.6g}, Bias={bias:.6g}")
+ predicts = np.ceil(intercept * ns[:, None] + bias - np.real(x))
+ print(f"Error: {np.mean(np.abs(bests - predicts))}")
diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb
deleted file mode 100644
index 44f3abd..0000000
--- a/buch/papers/laguerre/scripts/gamma_approx.ipynb
+++ /dev/null
@@ -1,395 +0,0 @@
-{
- "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/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
new file mode 100644
index 0000000..5b09e59
--- /dev/null
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -0,0 +1,116 @@
+from pathlib import Path
+
+import numpy as np
+import scipy.special
+
+EPSILON = 1e-7
+root = str(Path(__file__).parent)
+img_path = f"{root}/../images"
+fontsize = "medium"
+cmap = "plasma"
+
+
+def _prep_zeros_and_weights(x, w, n):
+ if x is None or w is None:
+ return np.polynomial.laguerre.laggauss(n)
+ return x, w
+
+
+def drop_imag(z):
+ if abs(z.imag) <= EPSILON:
+ z = z.real
+ return z
+
+
+def pochhammer(z, n):
+ return np.prod(z + np.arange(n))
+
+
+def find_optimal_shift(z, n):
+ mhat = 1.34093 * n + 0.854093
+ steps = int(np.floor(mhat - np.real(z)))
+ return steps
+
+
+def get_shifting_factor(z, steps):
+ factor = 1.0
+ if steps > 0:
+ factor = 1 / pochhammer(z, steps)
+ elif steps < 0:
+ factor = pochhammer(z + steps, -steps)
+ return factor
+
+
+def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
+ x, w = _prep_zeros_and_weights(x, w, n)
+ n = len(x)
+
+ z += 0j
+ steps = int(np.floor(target - np.real(z)))
+ z_shifted = z + steps
+ correction_factor = get_shifting_factor(z, steps)
+
+ res = np.sum(x ** (z_shifted - 1) * w)
+ res *= correction_factor
+ res = drop_imag(res)
+ return res
+
+
+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)
+
+ z += 0j
+ opt_shift = find_optimal_shift(z, n)
+ correction_factor = get_shifting_factor(z, opt_shift)
+ z_shifted = z + opt_shift
+
+ res = np.sum(x ** (z_shifted - 1) * w)
+ res *= correction_factor
+ res = drop_imag(res)
+ return res
+
+
+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)
+ res = drop_imag(res)
+ return res
+
+
+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:
+ return np.pi / (
+ np.sin(np.pi * z) * laguerre_gamma_simple(1 - z, x, w)
+ ) # Reflection formula
+ 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)
+ if func == "simple":
+ f = laguerre_gamma_simple
+ elif func == "mirror":
+ f = laguerre_gamma_mirror
+ elif func == "optimal_shifted":
+ f = laguerre_gamma_opt_shifted
+ else:
+ f = laguerre_gamma_shifted
+ return np.array([f(zi, x, w, n, **kwargs) for zi in z])
+
+
+def calc_rel_error(x, y):
+ 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/integrand.py b/buch/papers/laguerre/scripts/integrand.py
new file mode 100644
index 0000000..e970721
--- /dev/null
+++ b/buch/papers/laguerre/scripts/integrand.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""Plot for integrand of gamma function with shifting terms."""
+
+if __name__ == "__main__":
+ import os
+ from pathlib import Path
+
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ EPSILON = 1e-12
+ xlims = np.array([-3, 3])
+
+ root = str(Path(__file__).parent)
+ img_path = f"{root}/../images"
+ os.makedirs(img_path, exist_ok=True)
+
+ t = np.logspace(*xlims, 1001)[:, None]
+
+ z = np.array([-4.5, -2, -1, -0.5, 0.0, 0.5, 1, 2, 4.5])
+ r = t ** z
+
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(4, 2.4))
+ ax.semilogx(t, r)
+ ax.set_xlim(*(10.0 ** xlims))
+ ax.set_ylim(1e-3, 40)
+ ax.set_xlabel(r"$x$")
+ # ax.set_ylabel(r"$x^z$")
+ ax.grid(1, "both")
+ labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)]
+ ax.legend(labels, ncol=2, loc="upper left", fontsize="small")
+ fig.savefig(f"{img_path}/integrand.pdf")
diff --git a/buch/papers/laguerre/scripts/integrand_exp.py b/buch/papers/laguerre/scripts/integrand_exp.py
new file mode 100644
index 0000000..e649b26
--- /dev/null
+++ b/buch/papers/laguerre/scripts/integrand_exp.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""Plot for integrand of gamma function with shifting terms."""
+
+if __name__ == "__main__":
+ import os
+ from pathlib import Path
+
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ EPSILON = 1e-12
+ xlims = np.array([-3, 3])
+
+ root = str(Path(__file__).parent)
+ img_path = f"{root}/../images"
+ os.makedirs(img_path, exist_ok=True)
+
+ t = np.logspace(*xlims, 1001)[:, None]
+
+ z = np.array([-1, -0.5, 0.0, 0.5, 1, 2, 3, 4, 4.5])
+ e = np.exp(-t)
+ r = t ** z * e
+
+ fig, ax = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(4, 2.4))
+ ax.semilogx(t, r)
+ # ax.plot(t,np.exp(-t))
+ ax.set_xlim(10 ** (-2), 20)
+ ax.set_ylim(1e-3, 10)
+ ax.set_xlabel(r"$x$")
+ # ax.set_ylabel(r"$x^z e^{-x}$")
+ ax.grid(1, "both")
+ labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)]
+ ax.legend(labels, ncol=2, loc="upper left", fontsize="small")
+ # fig.savefig(f"{img_path}/integrand_exp.pgf")
+ fig.savefig(f"{img_path}/integrand_exp.pdf")
+ # plt.show()
diff --git a/buch/papers/laguerre/scripts/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py
deleted file mode 100644
index b9088d0..0000000
--- a/buch/papers/laguerre/scripts/laguerre_plot.py
+++ /dev/null
@@ -1,100 +0,0 @@
-#!/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/scripts/laguerre_poly.py b/buch/papers/laguerre/scripts/laguerre_poly.py
new file mode 100644
index 0000000..9700ab4
--- /dev/null
+++ b/buch/papers/laguerre/scripts/laguerre_poly.py
@@ -0,0 +1,108 @@
+import numpy as np
+
+
+def get_ticks(start, end, step=1):
+ ticks = np.arange(start, end, step)
+ return ticks[ticks != 0]
+
+
+if __name__ == "__main__":
+ import os
+ from pathlib import Path
+
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import scipy.special as ss
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ 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(get_ticks(-ylim, ylim), minor=True)
+ ax.set_yticks(get_ticks(-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_poly.pgf")
+ fig.savefig(f"{img_path}/laguerre_poly.pdf")
+ # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_complex.py b/buch/papers/laguerre/scripts/rel_error_complex.py
new file mode 100644
index 0000000..4a714fa
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_complex.py
@@ -0,0 +1,43 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+ import scipy.special
+
+ import gamma_approx as ga
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ xmax = 4
+ vals = np.linspace(-xmax + ga.EPSILON, xmax, 100)
+ x, y = np.meshgrid(vals, vals)
+ mesh = x + 1j * y
+ input = mesh.flatten()
+
+ lanczos = scipy.special.gamma(mesh)
+ lag = ga.eval_laguerre_gamma(input, n=8, func="optimal_shifted").reshape(mesh.shape)
+ rel_error = np.abs(ga.calc_rel_error(lanczos, lag))
+
+ fig, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(3.5, 2.1))
+ _c = ax.pcolormesh(
+ x, y, rel_error, shading="gouraud", cmap=ga.cmap, norm=mpl.colors.LogNorm()
+ )
+ cbr = fig.colorbar(_c, ax=ax)
+ cbr.minorticks_off()
+ # ax.set_title("Relative Error")
+ ax.set_xlabel("Re($z$)")
+ ax.set_ylabel("Im($z$)")
+ minor_ticks = np.arange(-xmax, xmax + ga.EPSILON)
+ ticks = np.arange(-xmax, xmax + ga.EPSILON, 2)
+ ax.set_xticks(ticks)
+ ax.set_xticks(minor_ticks, minor=True)
+ ax.set_yticks(ticks)
+ ax.set_yticks(minor_ticks, minor=True)
+ fig.savefig(f"{ga.img_path}/rel_error_complex.pdf")
+ # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_mirror.py b/buch/papers/laguerre/scripts/rel_error_mirror.py
new file mode 100644
index 0000000..7348d5e
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_mirror.py
@@ -0,0 +1,38 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+ import scipy.special
+
+ import gamma_approx as ga
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ xmin = -15
+ xmax = 15
+ ns = np.arange(2, 12, 2)
+ ylim = np.array([-11, 1])
+ x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400)
+ gamma = scipy.special.gamma(x)
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
+ for n in ns:
+ gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="mirror")
+ rel_err = ga.calc_rel_error(gamma, gamma_lag)
+ ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+ ax.set_xlim(x[0], x[-1])
+ ax.set_ylim(*(10.0 ** ylim))
+ ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5))
+ ax.set_xticks(np.arange(xmin, xmax), minor=True)
+ ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+ ax.set_xlabel(r"$z$")
+ # ax.set_ylabel("Relativer Fehler")
+ ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize)
+ ax.grid(1, "both")
+ # fig.savefig(f"{ga.img_path}/rel_error_mirror.pgf")
+ fig.savefig(f"{ga.img_path}/rel_error_mirror.pdf")
diff --git a/buch/papers/laguerre/scripts/rel_error_range.py b/buch/papers/laguerre/scripts/rel_error_range.py
new file mode 100644
index 0000000..ece3b6d
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_range.py
@@ -0,0 +1,41 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+ import scipy.special
+
+ import gamma_approx as ga
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+ N = 1201
+ xmax = 6
+ xmin = -xmax
+ ns = np.arange(2, 12, 2)
+ ylim = np.array([-11, -1.2])
+
+ x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, N)
+ gamma = scipy.special.gamma(x)
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2))
+ for n in ns:
+ gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+ rel_err = ga.calc_rel_error(gamma, gamma_lag)
+ ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+ ax.set_xlim(x[0], x[-1])
+ ax.set_ylim(*(10.0 ** ylim))
+ ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 2))
+ ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON), minor=True)
+ ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+ ax.set_yticks(10.0 ** np.arange(*ylim, 1), "", minor=True)
+ ax.set_xlabel(r"$z$")
+ # ax.set_ylabel("Relativer Fehler")
+ ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize)
+ ax.grid(1, "both")
+ # fig.savefig(f"{ga.img_path}/rel_error_range.pgf")
+ fig.savefig(f"{ga.img_path}/rel_error_range.pdf")
+ # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_shifted.py b/buch/papers/laguerre/scripts/rel_error_shifted.py
new file mode 100644
index 0000000..f53c89b
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_shifted.py
@@ -0,0 +1,40 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+ import scipy.special
+
+ import gamma_approx as ga
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+ n = 8 # order of Laguerre polynomial
+ N = 200 # number of points in interval
+
+ step = 1 / (N - 1)
+ x = np.linspace(step, 1 - step, N + 1)
+ targets = np.arange(10, 14)
+ gamma = scipy.special.gamma(x)
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.1))
+ for target in targets:
+ gamma_lag = ga.eval_laguerre_gamma(x, target=target, n=n, func="shifted")
+ rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lag))
+ ax.semilogy(x, rel_error, label=f"$m={target}$", linewidth=3)
+ gamma_lgo = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+ rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lgo))
+ ax.semilogy(x, rel_error, "m", linestyle=":", label="$m^*$", linewidth=3)
+ ax.set_xlim(x[0], x[-1])
+ ax.set_ylim(5e-9, 5e-8)
+ ax.set_xlabel(r"$z$")
+ ax.set_xticks(np.linspace(0, 1, 6))
+ ax.set_xticks(np.linspace(0, 1, 11), minor=True)
+ ax.grid(1, "both")
+ ax.legend(ncol=1, fontsize=ga.fontsize)
+ # fig.savefig(f"{ga.img_path}/rel_error_shifted.pgf")
+ fig.savefig(f"{ga.img_path}/rel_error_shifted.pdf")
+ # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_simple.py b/buch/papers/laguerre/scripts/rel_error_simple.py
new file mode 100644
index 0000000..686500b
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_simple.py
@@ -0,0 +1,41 @@
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+ import numpy as np
+ import scipy.special
+
+ import gamma_approx as ga
+
+ # mpl.rc("text", usetex=True)
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+ # mpl.rcParams.update({"font.family": "serif", "font.serif": "TeX Gyre Termes"})
+
+ # Simple / naive
+ xmin = -5
+ xmax = 30
+ ns = np.arange(2, 12, 2)
+ ylim = np.array([-11, 6])
+ x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400)
+ gamma = scipy.special.gamma(x)
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
+ for n in ns:
+ gamma_lag = ga.eval_laguerre_gamma(x, n=n)
+ rel_err = ga.calc_rel_error(gamma, gamma_lag)
+ ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+ ax.set_xlim(x[0], x[-1])
+ ax.set_ylim(*(10.0 ** ylim))
+ ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5))
+ ax.set_xticks(np.arange(xmin, xmax), minor=True)
+ ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+ ax.set_xlabel(r"$z$")
+ # ax.set_ylabel("Relativer Fehler")
+ ax.legend(ncol=3, fontsize=ga.fontsize)
+ ax.grid(1, "both")
+ # fig.savefig(f"{ga.img_path}/rel_error_simple.pgf")
+ fig.savefig(f"{ga.img_path}/rel_error_simple.pdf")
diff --git a/buch/papers/laguerre/scripts/targets.py b/buch/papers/laguerre/scripts/targets.py
new file mode 100644
index 0000000..3bc7f52
--- /dev/null
+++ b/buch/papers/laguerre/scripts/targets.py
@@ -0,0 +1,58 @@
+import numpy as np
+import scipy.special
+
+import gamma_approx as ga
+
+
+def find_best_loc(N=200, a=1.375, b=0.5, ns=None):
+ if ns is None:
+ ns = np.arange(2, 13)
+ bests = []
+ step = 1 / (N - 1)
+ x = np.linspace(step, 1 - step, N + 1)
+ gamma = scipy.special.gamma(x)
+ for n in ns:
+ zeros, weights = np.polynomial.laguerre.laggauss(n)
+ est = np.ceil(b + a * n)
+ targets = np.arange(max(est - 2, 0), est + 3)
+ rel_error = []
+ for target in targets:
+ gamma_lag = ga.eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted")
+ rel_error.append(np.abs(ga.calc_rel_error(gamma, gamma_lag)))
+ rel_error = np.stack(rel_error, -1)
+ best = np.argmin(rel_error, -1) + targets[0]
+ bests.append(best)
+ return np.stack(bests, 0)
+
+
+if __name__ == "__main__":
+ import matplotlib as mpl
+ import matplotlib.pyplot as plt
+
+ mpl.rcParams.update(
+ {
+ "mathtext.fontset": "stix",
+ "font.family": "serif",
+ "font.serif": "TeX Gyre Termes",
+ }
+ )
+
+ N = 200
+ ns = np.arange(2, 13)
+
+ bests = find_best_loc(N, ns=ns)
+
+ fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(3.5, 2.1))
+ v = ax.imshow(bests, cmap=ga.cmap, aspect="auto", interpolation="nearest")
+ plt.colorbar(v, ax=ax, label=r"$m^*$")
+ ticks = np.arange(0, N + 1, N // 5)
+ ax.set_xlim(0, 1)
+ ax.set_xticks(ticks)
+ ax.set_xticklabels([f"{v:.2f}" for v in ticks / N])
+ ax.set_xticks(np.arange(0, N + 1, N // 20), minor=True)
+ ax.set_yticks(np.arange(len(ns)))
+ ax.set_yticklabels(ns)
+ ax.set_xlabel(r"$z$")
+ ax.set_ylabel(r"$n$")
+ # fig.savefig(f"{ga.img_path}/targets.pgf")
+ fig.savefig(f"{ga.img_path}/targets.pdf")