From e1f5d6267540ea8dc758696fb08cb7540362cf8f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Mon, 18 Jul 2022 17:34:37 +0200
Subject: First complete draft of Laguerre chapter

---
 buch/papers/laguerre/gamma.tex | 242 +++++++++++++++++++++++++++--------------
 1 file changed, 163 insertions(+), 79 deletions(-)

(limited to 'buch/papers/laguerre/gamma.tex')

diff --git a/buch/papers/laguerre/gamma.tex b/buch/papers/laguerre/gamma.tex
index a28c180..eb64fa2 100644
--- a/buch/papers/laguerre/gamma.tex
+++ b/buch/papers/laguerre/gamma.tex
@@ -23,8 +23,8 @@ Integral der Form
 ,
 \quad
 \text{wobei Realteil von $z$ grösser als $0$}
-,
 \label{laguerre:gamma}
+.
 \end{align}
 Der Term $e^{-t}$ ist genau die Gewichtsfunktion der Laguerre-Integration und
 der Definitionsbereich passt ebenfalls genau für dieses Verfahren.
@@ -72,7 +72,7 @@ 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{Cassity1965AbcissasCA}
+dass die Berechnung der Gewichte $A_i$ nach \cite{laguerre:Cassity1965AbcissasCA}
 \begin{align*}
 A_i
 =
@@ -85,7 +85,7 @@ A_i
 }
 \end{align*}
 Evaluationen der Gamma-Funktion benötigen.
-Somit scheint diese Methode nicht geeignet für unser Vorhaben.
+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$.
@@ -95,10 +95,10 @@ 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 auf Grund der vorherigen Punkte,
+Somit entscheiden wir uns aufgrund der vorherigen Punkte,
 die zweite Variante weiterzuverfolgen.
 
-\subsubsection{Naiver Ansatz}
+\subsubsection{Direkter Ansatz}
 Wenden wir also die Gauss-Laguerre-Quadratur aus
 \eqref{laguerre:laguerrequadratur} auf die Gamma-Funktion
 \eqref{laguerre:gamma} an ergibt sich
@@ -111,15 +111,16 @@ Wenden wir also die Gauss-Laguerre-Quadratur aus
 
 \begin{figure}
 \centering
-\input{papers/laguerre/images/rel_error_simple.pgf}
-\vspace{-12pt}
-\caption{Relativer Fehler des naiven Ansatzes
+% \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}
 
 Bevor wir die Gauss-Laguerre-Quadratur anwenden,
-möchten wir als erstes eine Fehlerabschätzung durchführen.
+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
@@ -130,6 +131,7 @@ Für das Integral der Gamma-Funktion ergibt sich also
 \\
  & =
 (z - 2n)_{2n} \xi^{z - 2n - 1}
+.
 \end{align*}
 Eingesetzt im Fehlerterm \eqref{laguerre:lag_error} resultiert
 \begin{align}
@@ -147,17 +149,19 @@ 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.
 
-Wenden wir nun also naiv die Gauss-Laguerre-Quadratur auf die Gammafunktion an.
+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 < 2n$,
+wie der relative Fehler Nullstellen aufweist für ganzzahlige $z \leq 2n$,
 was laut der Theorie der Gauss-Quadratur auch zu erwarten ist,
 denn die Approximation via Gauss-Quadratur
-ist exakt für zu integrierende Polynome mit Grad $< 2n-1$.
+ist exakt für zu integrierende Polynome mit Grad $\leq 2n-1$
+und von $z$ auch 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.
@@ -168,9 +172,10 @@ könnten wir die Reflektionsformel der Gamma-Funktion ausnutzen.
 
 \begin{figure}
 \centering
-\input{papers/laguerre/images/rel_error_mirror.pgf}
-\vspace{-12pt}
-\caption{Relativer Fehler des naiven Ansatz mit Spiegelung negativer Realwerte
+% \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}
@@ -202,8 +207,9 @@ dadurch geeignete Gegenmassnahmen zu entwickeln.
 % und Abbildung~\ref{laguerre:fig:integrand_exp} grafisch dargestellt werden.
 \begin{figure}
 \centering
-\input{papers/laguerre/images/integrand.pgf}
-\vspace{-12pt}
+% \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}
@@ -211,7 +217,7 @@ dadurch geeignete Gegenmassnahmen zu entwickeln.
 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-
+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.
@@ -223,8 +229,9 @@ dass kleine Exponenten um $0$ genauere Resultate liefern sollten.
 
 \begin{figure}
 \centering
-\input{papers/laguerre/images/integrand_exp.pgf}
-\vspace{-12pt}
+% \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}
@@ -246,9 +253,9 @@ 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$ ist.
+grösser als $0$ und kleiner als $2n-1$ ist.
 
-\subsubsection{Finden der optimalen Berechnungsstelle}
+\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,
@@ -287,12 +294,13 @@ s(z, m)
 =
 \begin{cases}
 \displaystyle
-\frac{1}{(z - m)_m} & \text{wenn } m \geq 0 \\
-(z + m)_{-m}        & \text{wenn } m < 0
+\frac{1}{(z)_m} & \text{wenn } m \geq 0 \\
+(z + m)_{-m}    & \text{wenn } m < 0
 \end{cases}
 .
 \end{align*}
 
+\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,
@@ -305,9 +313,17 @@ 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
@@ -329,21 +345,14 @@ m^*
 \operatorname*{argmin}_m \max_\xi R_{n,m}(\xi)
 .
 \end{align*}
-Allerdings ist die Funktion $R_{n,m}(\xi)$ unbeschränkt.
+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 Sinn macht für kleine $n$,
+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.
 
-\begin{figure}
-\centering
-% \includegraphics{papers/laguerre/images/targets.pdf}
-\input{papers/laguerre/images/targets.pgf}
-\vspace{-12pt}
-\caption{$a$ in Abhängigkeit von $z$ und $n$}
-\label{laguerre:fig:targets}
-\end{figure}
-
 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,
@@ -369,11 +378,20 @@ Den linearen Regressor
 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}
+
 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 optimalen Verschiebungsterm
+Der optimalen Verschiebungsterm kann nun mit
 \begin{align*}
 m^*
 \approx
@@ -381,61 +399,127 @@ m^*
 =
 \lceil \alpha n + \beta - z \rceil
 \end{align*}
-kann nun mit dem linearen Regressor und $z$ gefunden werden.
-
-\begin{figure}
-\centering
-\input{papers/laguerre/images/estimates.pgf}
-\vspace{-12pt}
-\caption{Schätzung Mittelwert von $m$ und Fehler}
-\label{laguerre:fig:schaetzung}
-\end{figure}
-
-\subsection{Resultate}
-
-\subsubsection{Relativer Fehler}
+% 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}
-\vspace{-12pt}
+% \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,
+bei 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}
-\vspace{-12pt}
+% \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 Grade $n$ der Laguerre-Polynome}
+für verschiedene reele Werte von $z$ und Laguerre-Polynome vom Grad $n$}
 \label{laguerre:fig:rel_error_range}
 \end{figure}
 
-\subsubsection{Vergleich mit Lanczos-Methode}
-{\color{red}
-$ $\newline
-$n = 7$:\newline
-Lanczos Polynomgrad auf 13 Stellen.\newline
-Unsere Methode auf 7 Stellen
-}
-
-% 2. Die Fehlerabschätzung ist problematisch, 
-% weil die Funktion R_n(\xi) unbeschränkt ist.
-% Daher kann man nicht einfach nach dem Maximum von R_n(\xi) suchen.
-% Man muss zunächst irgendwie das \xi unter Kontrolle bringen. 
-% Das scheint mir äusserst schwierig zu sein.
+\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}
 
-% Ich möchte daher folgendes anregen: 
-% Im Sinne der Formulierung des Problems,
-% wie im Punkt 1 oben könnten Sie für verschiedene n
-% nach den optimalen Intervallen [a(n),a(n)+1] suchen,
-% und versuchen, einen empirischen Zusammenhang (Faustregel) 
-% zwischen n und a(n) zu formulieren.
-% Das ist etwa gleich gut,
-% da ja der Witz der Gauss-Integration ist,
-% dass man eben nur sehr kleine n überhaupt in Betracht zieht,
-% d.h. man braucht keine exakte Gesetzmässigkeit für a(n).
+\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$-$7$ 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$ Funktionasevaluationen, 
+wenigen Multiplikationen und Additionen besteht.
+Also könnte diese Methode z.B. Anwendung in Systemen mit wenig Rechenleistung
+und/oder knappen Energieressourcen finden.
\ No newline at end of file
-- 
cgit v1.2.1

