path: root/buch/papers/laguerre/gamma.tex

%
% 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 $\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{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}
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}
\vspace{-12pt}
\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}
\vspace{-12pt}
\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(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}
\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}
\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$ 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(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*}
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}
.
\end{align*}

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}
% 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.
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$,
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,
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.

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
\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.

\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}

\begin{figure}
\centering
\input{papers/laguerre/images/rel_error_shifted.pgf}
\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}
  
\begin{figure}
\centering
\input{papers/laguerre/images/rel_error_range.pgf}
\vspace{-12pt}
\caption{Relativer Fehler des Ansatzes mit optimalen Verschiebungsterm
für verschiedene reele Werte von $z$ und Grade $n$ der Laguerre-Polynome}
\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.

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