diff options
Diffstat (limited to 'buch')
40 files changed, 1912 insertions, 600 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..42cd6f6 100644 --- a/buch/papers/laguerre/definition.tex +++ b/buch/papers/laguerre/definition.tex @@ -18,11 +18,14 @@ x \in \mathbb{R} . \label{laguerre:dgl} \end{align} +Spannenderweise wurde die verallgemeinerte Laguerre-Differentialgleichung +zuerst von Yacovlevich Sonine (1849 - 1915) beschrieben, +aber auf Grund ihrer Ähnlichkeit wurde sie 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 +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. @@ -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 @@ -142,16 +154,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..9b901ae 100644 --- a/buch/papers/laguerre/eigenschaften.tex +++ b/buch/papers/laguerre/eigenschaften.tex @@ -3,35 +3,44 @@ % % (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. -} +% \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} +% Die Laguerre-Polynome besitzen einige interessante Eigenschaften +% \rhead{Eigenschaften} -\subsection{Orthogonalität - \label{laguerre:subsection:orthogonal}} +% \subsection{Orthogonalität +% \label{laguerre:subsection:orthogonal}} +\section{Orthogonalität + \label{laguerre:section:orthogonal}} Im Abschnitt~\ref{laguerre:section:definition} haben wir behauptet, dass die Laguerre-Polynome orthogonale Polynome sind. Zu dieser Behauptung möchten wir nun einen Beweis liefern. -Wenn wir die Laguerre\--Differentialgleichung in ein -Sturm\--Liouville\--Problem umwandeln können, haben wir bewiesen, dass es sich -bei -den Laguerre\--Polynomen um orthogonale Polynome handelt (siehe +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, +in dem wir diese Operatoren einander gleichsetzen. Aus der Beziehung \begin{align} S @@ -56,17 +65,19 @@ 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 +117,14 @@ 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} +Damit können wir schlussfolgern, dass die verallgemeinerten Laguerre-Polynome +orthogonal bezüglich des Skalarproduktes auf dem Intervall $(0, \infty)$ +mit der verallgemeinerten Laguerre\--Gewichtsfunktion $w(x)=x^\nu e^{-x}$ sind. +Die Laguerre-Polynome ($\nu=0$) sind somit orthognal im Intervall $(0, \infty)$ +mit der Gewichtsfunktion $w(x)=e^{-x}$. -\subsection{Drei-Terme Rekursion} +% \subsection{Rodrigues-Formel} -\subsection{Beziehung mit der Hypergeometrischen Funktion} +% \subsection{Drei-Terme Rekursion} +% \subsection{Beziehung mit der Hypergeometrischen Funktion} diff --git a/buch/papers/laguerre/gamma.tex b/buch/papers/laguerre/gamma.tex index e3838b0..b76daeb 100644 --- a/buch/papers/laguerre/gamma.tex +++ b/buch/papers/laguerre/gamma.tex @@ -19,58 +19,507 @@ 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}$ 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 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$, +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 $\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. +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 optimalen 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, +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} +\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$-$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 diff --git a/buch/papers/laguerre/images/estimates.pdf b/buch/papers/laguerre/images/estimates.pdf Binary files differnew file mode 100644 index 0000000..bd995de --- /dev/null +++ b/buch/papers/laguerre/images/estimates.pdf diff --git a/buch/papers/laguerre/images/gammaplot.pdf b/buch/papers/laguerre/images/gammaplot.pdf Binary files differnew file mode 100644 index 0000000..7c195f2 --- /dev/null +++ b/buch/papers/laguerre/images/gammaplot.pdf diff --git a/buch/papers/laguerre/images/integrand.pdf b/buch/papers/laguerre/images/integrand.pdf Binary files differnew file mode 100644 index 0000000..76be412 --- /dev/null +++ b/buch/papers/laguerre/images/integrand.pdf diff --git a/buch/papers/laguerre/images/integrand_exp.pdf b/buch/papers/laguerre/images/integrand_exp.pdf Binary files differnew file mode 100644 index 0000000..5fe7a7b --- /dev/null +++ b/buch/papers/laguerre/images/integrand_exp.pdf diff --git a/buch/papers/laguerre/images/laguerre_poly.pdf b/buch/papers/laguerre/images/laguerre_poly.pdf Binary files differnew file mode 100644 index 0000000..21278f5 --- /dev/null +++ b/buch/papers/laguerre/images/laguerre_poly.pdf diff --git a/buch/papers/laguerre/images/laguerre_polynomes.pdf b/buch/papers/laguerre/images/laguerre_polynomes.pdf Binary files differdeleted file mode 100644 index 3976bc7..0000000 --- a/buch/papers/laguerre/images/laguerre_polynomes.pdf +++ /dev/null diff --git a/buch/papers/laguerre/images/rel_error_complex.pdf b/buch/papers/laguerre/images/rel_error_complex.pdf Binary files differnew file mode 100644 index 0000000..c7bb37a --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_complex.pdf diff --git a/buch/papers/laguerre/images/rel_error_mirror.pdf b/buch/papers/laguerre/images/rel_error_mirror.pdf Binary files differnew file mode 100644 index 0000000..8a27d41 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_mirror.pdf diff --git a/buch/papers/laguerre/images/rel_error_range.pdf b/buch/papers/laguerre/images/rel_error_range.pdf Binary files differnew file mode 100644 index 0000000..bb8a2d7 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_range.pdf diff --git a/buch/papers/laguerre/images/rel_error_shifted.pdf b/buch/papers/laguerre/images/rel_error_shifted.pdf Binary files differnew file mode 100644 index 0000000..b7e72dc --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_shifted.pdf diff --git a/buch/papers/laguerre/images/rel_error_simple.pdf b/buch/papers/laguerre/images/rel_error_simple.pdf Binary files differnew file mode 100644 index 0000000..3212e42 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_simple.pdf diff --git a/buch/papers/laguerre/images/targets.pdf b/buch/papers/laguerre/images/targets.pdf Binary files differnew file mode 100644 index 0000000..9514a6d --- /dev/null +++ b/buch/papers/laguerre/images/targets.pdf diff --git a/buch/papers/laguerre/main.tex b/buch/papers/laguerre/main.tex index 00e3b43..d69fbed 100644 --- a/buch/papers/laguerre/main.tex +++ b/buch/papers/laguerre/main.tex @@ -8,7 +8,21 @@ \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 Approximationsmethoden +für das Integral $\int_0^\infty \exp(-x) / x \, dx$ 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, +eine geeignete Approximation für die Gamma-Funktion zu finden +mittels Laguerre-Polynomen und der Gauss-Quadratur. + +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 Binary files differnew file mode 100644 index 0000000..3d00de3 --- /dev/null +++ b/buch/papers/laguerre/presentation/presentation.pdf 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..75858df 100644 --- a/buch/papers/laguerre/quadratur.tex +++ b/buch/papers/laguerre/quadratur.tex @@ -5,27 +5,62 @@ % \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 ausnützt. +Herleitungen und Analysen der Gauss-Quadratur können im +Abschnitt~\ref{buch:orthogonalitaet: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 gut funktioniert, +sollte es auch für andere Funktionen nicht schlecht funktionieren. +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 eine gewichtete 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 Polynome mit einer Funktion multiplizieren, +die schneller als jedes Polynom gegen $0$ geht, +damit das Integral immer noch konvergiert. +Die Laguerre-Polynome $L_n$ bieten 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} @@ -33,7 +68,7 @@ Gleichung~\eqref{laguerre:laguerrequadratur} lässt sich wiefolgt umformulieren: 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. +als Gewichte $A_i$ die Integrale $l_i(x)e^{-x}$ verwendet werden. Dabei sind \begin{align*} l_i(x_j) @@ -41,39 +76,111 @@ 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") |