%
% gamma.tex 
%
% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule
%
\section{Anwendung: Berechnung der Gamma-Funktion
  \label{laguerre:section:quad-gamma}}
Die Gauss-Laguerre-Quadratur kann nun verwendet werden,
um exponentiell abfallende Funktionen im Definitionsbereich $(0, \infty)$ zu
berechnen.
Dabei bietet sich z.B. die Gamma-Funkion bestens an, wie wir in den folgenden
Abschnitten sehen werden.

\subsection{Gamma-Funktion}
Die Gamma-Funktion ist eine Erweiterung der Fakultät auf die reale und komplexe
Zahlenmenge.
Die Definition~\ref{buch:rekursion:def:gamma} beschreibt die Gamma-Funktion als
Integral der Form
\begin{align}
\Gamma(z)
 & =
\int_0^\infty x^{z-1} e^{-x} \, dx
,
\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.
Zu erwähnen ist auch, dass für die verallgemeinerte Laguerre-Integration die
Gewichtsfunktion $t^\nu e^{-t}$ genau dem Integranden für $\nu=z-1$ entspricht.

\subsubsection{Funktionalgleichung}
Die Gamma-Funktion besitzt die gleiche Rekursionsbeziehung wie die Fakultät,
nämlich
\begin{align}
z \Gamma(z)
=
\Gamma(z+1)
.
\label{laguerre:gamma_funktional}
\end{align}

\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 $\operatorname{Re} 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{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 scheint diese Methode 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 auf Grund der vorherigen Punkte,
die zweite Variante weiterzuverfolgen.

\subsubsection{Naiver Ansatz}

\begin{figure}
\centering
\input{papers/laguerre/images/rel_error_simple.pgf}
\caption{Relativer Fehler des naiven 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.
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}
,
\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.

Wenden wir nun also naiv die Gauss-Laguerre-Quadratur auf die Gammafunktion 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$,
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$.
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}
\caption{Relativer Fehler des naiven Ansatz 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, a+1]$, $a \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/integrands.pgf}
\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/integrands_exp.pgf}
\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$,
welches das Intervall $[a,a+1]$ definiert,
in dem der relative Fehler minimal ist,
grösser als $0$ und abhängig von $n$ ist.

\subsubsection{Finden der optimalen Berechnungsstelle}
% 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 in einem geeigneten Intervall $[a, a+1]$, 
$a \in \mathbb{Z}$, 
evaluiert und dann mit der 
Funktionalgleichung \eqref{laguerre:gamma_funktional} zurückverschiebt.
Aus Gründen der Übersichtlichkeit möchten wir das Problem nur für reele Zahlen
formulieren.

Die optimale Stelle $a \leq z^* \leq a+1$,
mit $z^* \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^*$ ist.
Wir fügen einen Verschiebungsterm um $m \in \mathbb{N}$ Stellen ein, 
daraus folgt
\begin{align*}
R_{n,m}(\xi)
=
\frac{(z - 2n)_{2n}}{(z - m)_m} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1}
,\quad
\text{für }
\xi \in (0, \infty)
,
\end{align*}
wobei $z^* = z + m - 1$ ist.
Das Optimierungsproblem daraus lässt sich als
\begin{align*}
\operatorname*{argmin}_m \max_\xi R_{n,m}(\xi)
\end{align*}
formulieren.
Allerdings ist die Funktion $R_{n,m}(\xi)$ unbeschränkt.
Dazu müssten wir $\xi$ versuchen unter Kontroller zu bringen,
was ein äussersts schwieriges Unterfangen zu sein scheint.
Da die Gauss-Quadratur aber sowieso nur wirklich Sinn macht für kleine $n$,
können die Intervalle $[a(n), a(n)+1]$ empirisch gesucht werden.


\begin{align*}
m^*
=
\lceil \alpha n + \beta + \lfloor z \rfloor - z \rceil - \lfloor z \rfloor
\end{align*}

\begin{figure}
\centering
\includegraphics{papers/laguerre/images/targets.pdf}
\caption{$a$ in Abhängigkeit von $z$ und $n$}
\label{laguerre:fig:targets}
\end{figure}


\begin{figure}
\centering
\input{papers/laguerre/images/estimate.pgf}
\caption{Schätzung Mittelwert von $m$ und Fehler}
\label{laguerre:fig:schaetzung}
\end{figure}
% 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.

% 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).


{
\large \color{red}
TODO:
Geeignete Minimierung für Fehler finden, so dass sie mit den emprisich
bestimmen optimalen Punkten übereinstimmen.
}