diff options
author | LordMcFungus <mceagle117@gmail.com> | 2022-07-22 21:28:45 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2022-07-22 21:28:45 +0200 |
commit | 23f17598c1742c70f442b94044a20aa821022c5a (patch) | |
tree | a945540ee6a4e86b37df2f01e3a91584b4797c4f /buch/papers/laguerre | |
parent | Merge pull request #2 from AndreasFMueller/master (diff) | |
parent | Merge pull request #25 from JODBaer/master (diff) | |
download | SeminarSpezielleFunktionen-23f17598c1742c70f442b94044a20aa821022c5a.tar.gz SeminarSpezielleFunktionen-23f17598c1742c70f442b94044a20aa821022c5a.zip |
Merge pull request #3 from AndreasFMueller/master
update
Diffstat (limited to '')
41 files changed, 3323 insertions, 148 deletions
diff --git a/buch/papers/laguerre/Makefile b/buch/papers/laguerre/Makefile index 606d7e1..85a1b83 100644 --- a/buch/papers/laguerre/Makefile +++ b/buch/papers/laguerre/Makefile @@ -3,7 +3,42 @@ # # (c) 2020 Prof Dr Andreas Mueller # +IMGFOLDER := images +PRESFOLDER := presentation -images: - @echo "no images to be created in laguerre" +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 +.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 1eb5034..39b5d6f 100644 --- a/buch/papers/laguerre/Makefile.inc +++ b/buch/papers/laguerre/Makefile.inc @@ -9,8 +9,5 @@ dependencies-laguerre = \ papers/laguerre/references.bib \ papers/laguerre/definition.tex \ papers/laguerre/eigenschaften.tex \ - papers/laguerre/quadratur.tex \ - papers/laguerre/transformation.tex \ - papers/laguerre/wasserstoff.tex - - + papers/laguerre/quadratur.tex \ + papers/laguerre/gamma.tex diff --git a/buch/papers/laguerre/definition.tex b/buch/papers/laguerre/definition.tex index 5f6d8bd..4729a93 100644 --- a/buch/papers/laguerre/definition.tex +++ b/buch/papers/laguerre/definition.tex @@ -4,45 +4,157 @@ % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % \section{Definition -\label{laguerre:section:definition}} + \label{laguerre:section:definition}} \rhead{Definition} - +Die verallgemeinerte Laguerre-Differentialgleichung ist gegeben durch \begin{align} - x y''(x) + (1 - x) y'(x) + n y(x) - = - 0 - \label{laguerre:dgl} +x y''(x) + (\nu + 1 - x) y'(x) + n y(x) += +0 +, \quad +n \in \mathbb{N}_0 +, \quad +x \in \mathbb{R} +\label{laguerre:dgl} +. \end{align} - +Spannenderweise wurde die verallgemeinerte Laguerre-Differentialgleichung +zuerst von Yacovlevich Sonine (1849 - 1915) beschrieben, +aber aufgrund ihrer Ähnlichkeit nach Laguerre benannt. +Die klassische Laguerre-Diffentialgleichung erhält man, wenn $\nu = 0$. +Hier wird die verallgemeinerte Laguerre-Differentialgleichung verwendet, +weil die Lösung mit derselben Methode berechnet werden kann. +Zusätzlich erhält man aber die Lösung für den allgmeinen Fall. +Zur Lösung von \eqref{laguerre:dgl} verwenden wir einen +Potenzreihenansatz. +Da wir bereits wissen, dass die Lösung orthogonale Polynome sind, +erscheint dieser Ansatz sinnvoll. +Setzt man nun den Ansatz +\begin{align*} +y(x) + & = +\sum_{k=0}^\infty a_k x^k +\\ +y'(x) + & = +\sum_{k=1}^\infty k a_k x^{k-1} += +\sum_{k=0}^\infty (k+1) a_{k+1} x^k +\\ +y''(x) + & = +\sum_{k=2}^\infty k (k-1) a_k x^{k-2} += +\sum_{k=1}^\infty (k+1) k a_{k+1} x^{k-1} +\end{align*} +in die Differentialgleichung ein, erhält man +\begin{align*} +\sum_{k=1}^\infty (k+1) k a_{k+1} x^k ++ +(\nu + 1)\sum_{k=0}^\infty (k+1) a_{k+1} x^k +- +\sum_{k=0}^\infty k a_k x^k ++ +n \sum_{k=0}^\infty a_k x^k + & = +0 \\ +\sum_{k=1}^\infty +\left[ (k+1) k a_{k+1} + (\nu + 1)(k+1) a_{k+1} - k a_k + n a_k \right] x^k + & = +0. +\end{align*} +Daraus lässt sich die Rekursionsbeziehung +\begin{align*} +a_{k+1} + & = +\frac{k-n}{(k+1) (k + \nu + 1)} a_k +\end{align*} +ableiten. +Für ein konstantes $n$ erhalten wir als Potenzreihenlösung ein Polynom vom Grad +$n$, +denn für $k=n$ wird $a_{n+1} = 0$ und damit auch $a_{n+2}=a_{n+3}=\ldots=0$. +Aus der Rekursionsbeziehung ist zudem ersichtlich, +dass $a_0 \neq 0$ beliebig gewählt werden kann. +Wählen wir nun $a_0 = 1$, dann folgt für die Koeffizienten $a_1, a_2, a_3$ +\begin{align*} +a_1 += +-\frac{n}{1 \cdot (\nu + 1)} +, & & +a_2 += +\frac{(n-1)n}{1 \cdot 2 \cdot (\nu + 1)(\nu + 2)} +, & & +a_3 += +-\frac{(n-2)(n-1)n}{1 \cdot 2 \cdot 3 \cdot (\nu + 1)(\nu + 2)(\nu + 3)} +\end{align*} +und allgemein +\begin{align*} +k + & \leq +n: + & +a_k + & = +(-1)^k \frac{n!}{(n-k)!} \frac{1}{k!(\nu + 1)_k} += +\frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k} +\\ +k & >n: + & +a_k + & = +0. +\end{align*} +Somit erhalten wir für $\nu = 0$ die Laguerre-Polynome \begin{align} - L_n(x) - = - \sum_{k=0}^{n} - \frac{(-1)^k}{k!} - \begin{pmatrix} - n \\ - k - \end{pmatrix} - x^k - \label{laguerre:polynom} +L_n(x) += +\sum_{k=0}^{n} \frac{(-1)^k}{k!} \binom{n}{k} x^k +\label{laguerre:polynom} \end{align} - +und mit $\nu \in \mathbb{R}$ die verallgemeinerten Laguerre-Polynome \begin{align} - x y''(x) + (\alpha + 1 - x) y'(x) + n y(x) - = - 0 - \label{laguerre:generell_dgl} +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} -\begin{align} - L_n^\alpha (x) - = - \sum_{k=0}^{n} - \frac{(-1)^k}{k!} - \begin{pmatrix} - n + \alpha \\ - n - k - \end{pmatrix} - x^k - \label{laguerre:polynom} -\end{align} +\subsection{Analytische Fortsetzung} +Durch die analytische Fortsetzung erhalten wir zudem noch die zweite Lösung der +Differentialgleichung mit der Form +\begin{align*} +\Xi_n(x) += +L_n(x) \ln(x) + \sum_{k=1}^\infty d_k x^k +. +\end{align*} +Nach einigen aufwändigen Rechnungen, +% die am besten ein Computeralgebrasystem übernimmt, +die den Rahmen dieses Kapitel sprengen würden, +erhalten wir +\begin{align*} +\Xi_n += +L_n(x) \ln(x) ++ +\sum_{k=1}^n \frac{(-1)^k}{k!} \binom{n}{k} +(\alpha_{n-k} - \alpha_n - 2 \alpha_k)x^k ++ +(-1)^n \sum_{k=1}^\infty \frac{(k-1)!n!}{((n+k)!)^2} x^{n+k}, +\end{align*} +wobei $\alpha_0 = 0$ und $\alpha_k =\sum_{i=1}^k i^{-1}$, +$\forall k \in \mathbb{N}$. +% 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 b7597e5..4adbe86 100644 --- a/buch/papers/laguerre/eigenschaften.tex +++ b/buch/papers/laguerre/eigenschaften.tex @@ -3,6 +3,111 @@ % % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % -\section{Eigenschaften -\label{laguerre:section:eigenschaften}} -\rhead{Eigenschaften}
\ No newline at end of file +\section{Orthogonalität + \label{laguerre:section:orthogonal}} +Im Abschnitt~\ref{laguerre:section:definition} +haben wir die Behauptung aufgestellt, +dass die Laguerre-Polynome orthogonal sind. +Zu dieser Behauptung möchten wir nun einen Beweis liefern. +Wenn wir \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 Beweis kann äquivalent auch über den Sturm-Liouville-Operator +\begin{align} +S += +\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right). +\label{laguerre:slop} +\end{align} +und den Laguerre-Operator +\begin{align} +\Lambda += +x \frac{d}{dx^2} + (\nu + 1 -x) \frac{d}{dx} +\end{align} +erhalten werden, +indem wir diese Operatoren einander gleichsetzen. +Aus der Beziehung +\begin{align} +S + & = +\Lambda +\nonumber +\\ +\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right) + & = +x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx} +\label{laguerre:sl-lag} +\end{align} +lässt sich sofort erkennen, dass $q(x) = 0$. +Ausserdem ist ersichtlich, dass $p(x)$ die Differentialgleichung +\begin{align*} +x \frac{dp}{dx} += +-(\nu + 1 - x) p +\end{align*} +erfüllen muss. +Durch Separation erhalten wir dann +\begin{align*} +\int \frac{dp}{p} + & = +-\int \frac{\nu + 1 - x}{x} \, dx += +-\int \frac{\nu + 1}{x} \, dx - \int 1\, dx +\\ +\log p + & = +-(\nu + 1)\log x - x + c +\\ +p(x) + & = +-C x^{\nu + 1} e^{-x} +. +\end{align*} +Eingefügt in Gleichung~\eqref{laguerre:sl-lag} ergibt sich +\begin{align*} +\frac{C}{w(x)} +\left( +x^{\nu+1} e^{-x} \frac{d^2}{dx^2} + +(\nu + 1 - x) x^{\nu} e^{-x} \frac{d}{dx} +\right) += +x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx}. +\end{align*} +Mittels Koeffizientenvergleich kann nun abgelesen werden, dass $w(x) = x^\nu +e^{-x}$ und $C=1$ mit $\nu > -1$. +Die Gewichtsfunktion $w(x)$ wächst für $x\rightarrow-\infty$ sehr schnell an, +deshalb ist die Laguerre-Gewichtsfunktion nur geeignet für den +Definitionsbereich $(0, \infty)$. +Bleibt nur noch sicherzustellen, dass die Randbedingungen, +\begin{align} +k_0 y(0) + h_0 p(0)y'(0) + & = +0 +\label{laguerre:sllag_randa} +\\ +k_\infty y(\infty) + h_\infty p(\infty) y'(\infty) + & = +0 +\label{laguerre:sllag_randb} +\end{align} +mit $|k_i|^2 + |h_i|^2 \neq 0,\,\forall i \in \{0, \infty\}$, erfüllt sind. +Am linken Rand (Gleichung~\eqref{laguerre:sllag_randa}) kann $y(0) = 1$, $k_0 = +0$ und $h_0 = 1$ verwendet werden, +was auch die Laguerre-Polynome ergeben haben. +Für den rechten Rand ist die Bedingung (Gleichung~\eqref{laguerre:sllag_randb}) +\begin{align*} +\lim_{x \rightarrow \infty} p(x) y'(x) + & = +\lim_{x \rightarrow \infty} -x^{\nu + 1} e^{-x} y'(x) += +0 +\end{align*} +für beliebige Polynomlösungen erfüllt für $k_\infty=0$ und $h_\infty=1$. +Damit können wir schlussfolgern: +Die verallgemeinerten Laguerre-Polynome sind orthogonal +bezüglich des Skalarproduktes auf dem Intervall $(0, \infty)$ +mit der verallgemeinerten Laguerre\--Gewichtsfunktion $w(x)=x^\nu e^{-x}$. +Die Laguerre-Polynome ($\nu=0$) sind somit orthognal im Intervall $(0, \infty)$ +mit der Gewichtsfunktion $w(x)=e^{-x}$. diff --git a/buch/papers/laguerre/gamma.tex b/buch/papers/laguerre/gamma.tex new file mode 100644 index 0000000..2e5fc06 --- /dev/null +++ b/buch/papers/laguerre/gamma.tex @@ -0,0 +1,532 @@ +% +% 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 hervorragend an, +wie wir in den folgenden Abschnitten sehen werden. + +\subsection{Gamma-Funktion} +Die Gamma-Funktion ist eine Erweiterung der Fakultät auf die reale und komplexe +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}$ im Integranden und der Integrationsbereich erfüllen +genau die Bedingungen der Laguerre-Integration. +% Der Term $e^{-t}$ ist genau die Gewichtsfunktion der Laguerre-Integration und +% der Definitionsbereich passt ebenfalls genau für dieses Verfahren. +Weiter zu erwähnen ist, dass für die verallgemeinerte Laguerre-Integration die +Gewichtsfunktion $t^\nu e^{-t}$ exakt dem Integranden für $\nu=z-1$ entspricht. + +\subsubsection{Funktionalgleichung} +Die 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{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} + +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} +, +\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 direkt die Gauss-Laguerre-Quadratur auf die Gamma-Funktion +an. +Dazu benötigen wir die Gewichte nach +\eqref{laguerre:quadratur_gewichte} +und als Stützstellen die Nullstellen des Laguerre-Polynomes $L_n$. +Evaluieren wir den relativen Fehler unserer Approximation zeigt sich ein +Bild wie in Abbildung~\ref{laguerre:fig:rel_error_simple}. +Man kann sehen, +wie der relative Fehler Nullstellen aufweist für ganzzahlige $z \leq 2n$. +Laut der Theorie der Gauss-Quadratur auch ist das zu erwarten, +da die Approximation via Gauss-Quadratur +exakt ist für zu integrierende Polynome mit Grad $\leq 2n-1$ +und hinzukommt, +dass zudem von $z$ noch $1$ abgezogen wird im Exponenten. +Es ist ersichtlich, +dass sich für den Polynomgrad $n$ ein Interval gibt, +in dem der relative Fehler minimal ist. +Links steigt der relative Fehler besonders stark an, +während er auf der rechten Seite zu konvergieren scheint. +Um die linke Hälfte in den Griff zu bekommen, +könnten wir die Reflektionsformel der Gamma-Funktion ausnutzen. + +\begin{figure} +\centering +% \input{papers/laguerre/images/rel_error_mirror.pgf} +\includegraphics{papers/laguerre/images/rel_error_mirror.pdf} +%\vspace{-12pt} +\caption{Relativer Fehler des Ansatzes mit Spiegelung negativer Realwerte +für verschiedene reele Werte von $z$ und Grade $n$ der Laguerre-Polynome} +\label{laguerre:fig:rel_error_mirror} +\end{figure} + +Spiegelt man nun $z$ mit negativem Realteil mittels der Reflektionsformel, +ergibt sich ein stabilerer Fehler in der linken Hälfte, +wie in Abbildung~\ref{laguerre:fig:rel_error_mirror}. +Die Spiegelung bringt nur für wenige Werte einen, +für praktische Anwendungen geeigneten, +relativen Fehler. +Wie wir aber in Abbildung~\ref{laguerre:fig:rel_error_simple} sehen konnten, +gibt es für jeden Polynomgrad $n$ ein Intervall $[a(n), a(n) + 1]$, +$a(n) \in \mathbb{Z}$, +in welchem der relative Fehler minimal ist. +Die Funktionalgleichung der Gamma-Funktion \eqref{laguerre:gamma_funktional} +könnte uns hier helfen, +das Problem in den Griff zu bekommen. + +\subsubsection{Analyse des Integranden} +Wie wir im vorherigen Abschnitt gesehen haben, +scheint der Integrand problematisch. +Darum möchten wir jetzt den Integranden analysieren, +um ihn besser verstehen zu können und +dadurch geeignete Gegenmassnahmen zu entwickeln. + +% Dieser Abschnitt soll eine grafisches Verständnis dafür schaffen, +% wieso der Integrand so problematisch ist. +% Was das heisst sollte in Abbildung~\ref{laguerre:fig:integrand} +% und Abbildung~\ref{laguerre:fig:integrand_exp} grafisch dargestellt werden. +\begin{figure} +\centering +% \input{papers/laguerre/images/integrand.pgf} +\includegraphics{papers/laguerre/images/integrand.pdf} +%\vspace{-12pt} +\caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$} +\label{laguerre:fig:integrand} +\end{figure} + +In Abbildung~\ref{laguerre:fig:integrand} ist der Integrand $x^z$ für +unterschiedliche Werte von $z$ dargestellt. +Dies entspricht der zu integrierenden Funktion $f(x)$ +der Gauss-Laguerre-Quadratur für die Gamma-Funktion. +Man erkennt, +dass für kleine $z$ sich ein singulärer Integrand ergibt +und auch für grosse $z$ wächst der Integrand sehr schnell an. +Das heisst, +die Ableitungen im Fehlerterm divergieren noch schneller +und das wirkt sich negativ auf die Genauigkeit der Approximation aus. +Somit lässt sich hier sagen, +dass kleine Exponenten um $0$ genauere Resultate liefern sollten. + +\begin{figure} +\centering +% \input{papers/laguerre/images/integrand_exp.pgf} +\includegraphics{papers/laguerre/images/integrand_exp.pdf} +%\vspace{-12pt} +\caption{Integrand $x^z e^{-x}$ mit unterschiedlichen Werten für $z$} +\label{laguerre:fig:integrand_exp} +\end{figure} + +In Abbildung~\ref{laguerre:fig:integrand_exp} fügen wir +die Dämpfung der Gewichtsfunktion $w(x)$ +der Gauss-Laguerre-Quadratur wieder hinzu +und erhalten so wieder den kompletten Integranden $x^{z-1} e^{-x}$ +der Gamma-Funktion. +Für negative $z$ ergeben sich immer noch Singularitäten, +wenn $x \rightarrow 0$. +Um $1$ wächst der Term $x^z$ schneller als die Dämpfung $e^{-x}$, +aber für $x \rightarrow \infty$ geht der Integrand gegen $0$. +Das führt zu glockenförmigen Kurven, +die für grosse Exponenten $z$ nach der Stelle $x=1$ schnell anwachsen. +Zu grosse Exponenten $z$ sind also immer noch problematisch. +Kleine positive $z$ scheinen nun also auch zulässig zu sein. +Damit formulieren wir die Vermutung, +dass $a(n)$, +welches das Intervall $[a(n), a(n) + 1]$ definiert, +in dem der relative Fehler minimal ist, +grösser als $0$ und kleiner als $2n-1$ ist. + +\subsubsection{Ansatz mit Verschiebungsterm} +% Mittels der Funktionalgleichung \eqref{laguerre:gamma_funktional} +% kann der Wert von $\Gamma(z)$ im Interval $z \in [a,a+1]$, +% in dem der relative Fehler minimal ist, +% evaluiert werden und dann mit der Funktionalgleichung zurückverschoben werden. +Nun stellt sich die Frage, +ob die Approximation mittels Gauss-Laguerre-Quadratur verbessert werden kann, +wenn man das Problem 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} & \text{wenn } m \geq 0 \\ +(z + m)_{-m} & \text{wenn } m < 0 +\end{cases} +. +\end{align*} + +\subsubsection{Finden der optimalen Berechnungsstelle} +Um die optimale Stelle $z^*(n) \in \left[a(n), a(n) + 1\right]$, +$z^*(n) \in \mathbb{R}$, +zu finden, +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} + +In Abbildung~\ref{laguerre:fig:schaetzung} sind die Resultate +der linearen Regression aufgezeigt mit $\alpha = 1.34094$ und $\beta = +0.854093$. +Die lineare Beziehung ist ganz klar ersichtlich und der Fit scheint zu genügen. +Der optimale Verschiebungsterm kann nun mit +\begin{align*} +m^* +\approx +\lceil \hat{m} - z \rceil += +\lceil \alpha n + \beta - z \rceil +\end{align*} +% kann nun mit dem linearen Regressor und $z$ +gefunden werden. + +\subsubsection{Evaluation des Schätzers} +In einem ersten Schritt möchten wir analysieren, +wie gut die Abschätzung des optimalen Verschiebungsterms ist. +Dazu bestimmen wir den relativen Fehler für verschiedene Verschiebungsterme $m$ +rund um $m^*$ bei gegebenem Polynomgrad $n = 8$ für $z \in (0, 1)$. +Abbildung~\ref{laguerre:fig:rel_error_shifted} sind die relativen Fehler +der Approximation dargestellt. +Man kann deutlich sehen, +dass der relative Fehler anwächst, +je weiter der Verschiebungsterm vom idealen Wert abweicht. +Zudem scheint der Schätzer den optimalen Verschiebungsterm gut zu bestimmen, +da der Schätzer zuerst der grünen Linie folgt und +dann beim Übergang auf die orange Linie wechselt. +\begin{figure} +\centering +% \input{papers/laguerre/images/rel_error_shifted.pgf} +\includegraphics{papers/laguerre/images/rel_error_shifted.pdf} +%\vspace{-12pt} +\caption{Relativer Fehler des Ansatzes mit Verschiebungsterm +für verschiedene reele Werte von $z$ und Verschiebungsterme $m$. +Das verwendete Laguerre-Polynom besitzt den Grad $n = 8$. +$m^*$ bezeichnet hier den optimalen Verschiebungsterm.} +\label{laguerre:fig:rel_error_shifted} +\end{figure} + +\subsubsection{Resultate} +Das Verfahren scheint für den Grad $n=8$ und $z \in (0,1)$ gut zu funktioneren. +Es stellt sich nun die Frage, +wie der relative Fehler sich für verschiedene $z$ und $n$ verhält. +In Abbildung~\ref{laguerre:fig:rel_error_range} sind die relativen Fehler für +unterschiedliche $n$ dargestellt. +Der relative Fehler scheint immer noch Nullstellen aufzuweisen +für ganzzahlige $z$. +Durch das Verschieben ergibt sich jetzt aber, +wie zu erwarten war, +ein periodischer relativer Fehler mit einer Periodendauer von $1$. +Zudem lässt sich erkennen, +dass der Fehler abhängig von der Ordnung $n$ +des verwendeten Laguerre-Polynoms ist. +Wenn der Grad $n$ um $1$ erhöht wird, +verbessert sich die Genauigkeit des Resultats um etwa eine signifikante Stelle. + +In Abbildung~\ref{laguerre:fig:rel_error_complex} +ist der Betrag des relativen Fehlers in der komplexen Ebene dargestellt. +Je stärker der Imaginäranteil von $z$ von $0$ abweicht, +umso schlechter wird die Genauigkeit der Approximation. +Das erstaunt nicht weiter, +da die Gauss-Quadratur eigentlich nur für reelle Zahlen definiert ist. +Wenn der Imaginäranteil von $z$ ungefähr $0$ ist, +lässt sich das gleiche Bild beobachten wie in +Abbildung~\ref{laguerre:fig:rel_error_range}. + +\begin{figure} +\centering +% \input{papers/laguerre/images/rel_error_range.pgf} +\includegraphics{papers/laguerre/images/rel_error_range.pdf} +%\vspace{-12pt} +\caption{Relativer Fehler des Ansatzes mit optimalen Verschiebungsterm +für verschiedene reele Werte von $z$ und Laguerre-Polynome vom Grad $n$} +\label{laguerre:fig:rel_error_range} +\end{figure} + +\begin{figure} +\centering +\includegraphics{papers/laguerre/images/rel_error_complex.pdf} +%\vspace{-12pt} +\caption{Absolutwert des relativen Fehlers in der komplexen Ebene} +\label{laguerre:fig:rel_error_complex} +\end{figure} + +\subsubsection{Vergleich mit Lanczos-Methode} +Nun stellt sich die Frage, +wie das in diesem Abschnitt beschriebene Approximationsverfahren +der Gamma-Funktion sich gegenüber den üblichen Approximationsverfahren schlägt. +Eine häufig verwendete Methode ist die Lanczos-Approximation, +welche gegeben ist durch +\begin{align} +\Gamma(z + 1) +\approx +\sqrt{2\pi} \left( z + \sigma + \frac{1}{2} \right)^{z + 1/2} +e^{-(z + \sigma + 1/2)} \sum_{k=0}^n g_k H_k(z) +, +\end{align} +wobei +\begin{align*} +g_k = \frac{e^\sigma \varepsilon_k (-1)^k}{\sqrt{2\pi}} +\sum_{r=0}^k (-1)^r \, \binom{k}{r} \, (k)_r +\left( \frac{e}{r + \sigma + \frac{1}{2}}\right)^{r + 1/2} +, +\end{align*} +\begin{align*} +\varepsilon_k += +\begin{cases} +1 & \text{für } k = 0 \\ +2 & \text{sonst} +\end{cases} +\quad \text{und}\quad +H_k(z) += +\frac{(-1)^k (-z)_k}{(z+1)_k} +\end{align*} +mit $H_0 = 1$ und $\sum_0^n g_k = 1$ (siehe \cite{laguerre:lanczos}). +Diese Methode wurde zum Beispiel in +{\em GNU Scientific Library}, {\em Boost}, {\em CPython} und +{\em musl} implementiert. +Diese Methode erreicht für $n = 7$ typischerweise Genauigkeit von $13$ +korrekten, signifikanten Stellen für reele Argumente. +Zum Vergleich: die vorgestellte Methode erreicht für $n = 7$ +eine minimale Genauigkeit von $6$ korrekten, signifikanten Stellen +für reele Argumente. +Das Resultat ist etwas enttäuschend, +aber nicht unerwartet, +da die Lanczos-Methode spezifisch auf dieses Problem zugeschnitten ist und +unsere Methode eine erweiterte allgemeine Methode ist. +Was die Komplexität der Berechnungen im Betrieb angeht, +ist die Gauss-Laguerre-Quadratur wesentlich ressourcensparender, +weil sie nur aus $n$ Funktionsevaluationen, +wenigen Multiplikationen und Additionen besteht. +Demzufolge könnte diese Methode Anwendung in Systemen mit wenig Rechenleistung +und/oder knappen Energieressourcen finden.
\ No newline at end of file diff --git a/buch/papers/laguerre/images/estimates.pdf b/buch/papers/laguerre/images/estimates.pdf 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/gammapaths.tex b/buch/papers/laguerre/images/gammapaths.tex new file mode 100644 index 0000000..efa0863 --- /dev/null +++ b/buch/papers/laguerre/images/gammapaths.tex @@ -0,0 +1,1024 @@ +\def\gammaplus{({\dx*0.0190},{\dy*52.0728}) + -- ({\dx*0.0200},{\dy*49.4422}) + -- ({\dx*0.0400},{\dy*24.4610}) + -- ({\dx*0.0600},{\dy*16.1457}) + -- ({\dx*0.0800},{\dy*11.9966}) + -- ({\dx*0.1000},{\dy*9.5135}) + -- ({\dx*0.1200},{\dy*7.8633}) + -- ({\dx*0.1400},{\dy*6.6887}) + -- ({\dx*0.1600},{\dy*5.8113}) + -- ({\dx*0.1800},{\dy*5.1318}) + -- ({\dx*0.2000},{\dy*4.5908}) + -- ({\dx*0.2200},{\dy*4.1505}) + -- ({\dx*0.2400},{\dy*3.7855}) + -- ({\dx*0.2600},{\dy*3.4785}) + -- ({\dx*0.2800},{\dy*3.2169}) + -- ({\dx*0.3000},{\dy*2.9916}) + -- ({\dx*0.3200},{\dy*2.7958}) + -- ({\dx*0.3400},{\dy*2.6242}) + -- ({\dx*0.3600},{\dy*2.4727}) + -- ({\dx*0.3800},{\dy*2.3383}) + -- ({\dx*0.4000},{\dy*2.2182}) + -- ({\dx*0.4200},{\dy*2.1104}) + -- ({\dx*0.4400},{\dy*2.0132}) + -- ({\dx*0.4600},{\dy*1.9252}) + -- ({\dx*0.4800},{\dy*1.8453}) + -- ({\dx*0.5000},{\dy*1.7725}) + -- ({\dx*0.5200},{\dy*1.7058}) + -- ({\dx*0.5400},{\dy*1.6448}) + -- ({\dx*0.5600},{\dy*1.5886}) + -- ({\dx*0.5800},{\dy*1.5369}) + -- ({\dx*0.6000},{\dy*1.4892}) + -- ({\dx*0.6200},{\dy*1.4450}) + -- ({\dx*0.6400},{\dy*1.4041}) + -- ({\dx*0.6600},{\dy*1.3662}) + -- ({\dx*0.6800},{\dy*1.3309}) + -- ({\dx*0.7000},{\dy*1.2981}) + -- ({\dx*0.7200},{\dy*1.2675}) + -- ({\dx*0.7400},{\dy*1.2390}) + -- ({\dx*0.7600},{\dy*1.2123}) + -- ({\dx*0.7800},{\dy*1.1875}) + -- ({\dx*0.8000},{\dy*1.1642}) + -- ({\dx*0.8200},{\dy*1.1425}) + -- ({\dx*0.8400},{\dy*1.1222}) + -- ({\dx*0.8600},{\dy*1.1031}) + -- ({\dx*0.8800},{\dy*1.0853}) + -- ({\dx*0.9000},{\dy*1.0686}) + -- ({\dx*0.9200},{\dy*1.0530}) + -- ({\dx*0.9400},{\dy*1.0384}) + -- ({\dx*0.9600},{\dy*1.0247}) + -- ({\dx*0.9800},{\dy*1.0119}) + -- ({\dx*1.0000},{\dy*1.0000}) + -- ({\dx*1.0200},{\dy*0.9888}) + -- ({\dx*1.0400},{\dy*0.9784}) + -- ({\dx*1.0600},{\dy*0.9687}) + -- ({\dx*1.0800},{\dy*0.9597}) + -- ({\dx*1.1000},{\dy*0.9514}) + -- ({\dx*1.1200},{\dy*0.9436}) + -- ({\dx*1.1400},{\dy*0.9364}) + -- ({\dx*1.1600},{\dy*0.9298}) + -- ({\dx*1.1800},{\dy*0.9237}) + -- ({\dx*1.2000},{\dy*0.9182}) + -- ({\dx*1.2200},{\dy*0.9131}) + -- ({\dx*1.2400},{\dy*0.9085}) + -- ({\dx*1.2600},{\dy*0.9044}) + -- ({\dx*1.2800},{\dy*0.9007}) + -- ({\dx*1.3000},{\dy*0.8975}) + -- ({\dx*1.3200},{\dy*0.8946}) + -- ({\dx*1.3400},{\dy*0.8922}) + -- ({\dx*1.3600},{\dy*0.8902}) + -- ({\dx*1.3800},{\dy*0.8885}) + -- ({\dx*1.4000},{\dy*0.8873}) + -- ({\dx*1.4200},{\dy*0.8864}) + -- ({\dx*1.4400},{\dy*0.8858}) + -- ({\dx*1.4600},{\dy*0.8856}) + -- ({\dx*1.4800},{\dy*0.8857}) + -- ({\dx*1.5000},{\dy*0.8862}) + -- ({\dx*1.5200},{\dy*0.8870}) + -- ({\dx*1.5400},{\dy*0.8882}) + -- ({\dx*1.5600},{\dy*0.8896}) + -- ({\dx*1.5800},{\dy*0.8914}) + -- ({\dx*1.6000},{\dy*0.8935}) + -- ({\dx*1.6200},{\dy*0.8959}) + -- ({\dx*1.6400},{\dy*0.8986}) + -- ({\dx*1.6600},{\dy*0.9017}) + -- ({\dx*1.6800},{\dy*0.9050}) + -- ({\dx*1.7000},{\dy*0.9086}) + -- ({\dx*1.7200},{\dy*0.9126}) + -- ({\dx*1.7400},{\dy*0.9168}) + -- ({\dx*1.7600},{\dy*0.9214}) + -- ({\dx*1.7800},{\dy*0.9262}) + -- ({\dx*1.8000},{\dy*0.9314}) + -- ({\dx*1.8200},{\dy*0.9368}) + -- ({\dx*1.8400},{\dy*0.9426}) + -- ({\dx*1.8600},{\dy*0.9487}) + -- ({\dx*1.8800},{\dy*0.9551}) + -- ({\dx*1.9000},{\dy*0.9618}) + -- ({\dx*1.9200},{\dy*0.9688}) + -- ({\dx*1.9400},{\dy*0.9761}) + -- ({\dx*1.9600},{\dy*0.9837}) + -- ({\dx*1.9800},{\dy*0.9917}) + -- ({\dx*2.0000},{\dy*1.0000}) + -- ({\dx*2.0200},{\dy*1.0086}) + -- ({\dx*2.0400},{\dy*1.0176}) + -- ({\dx*2.0600},{\dy*1.0269}) + -- ({\dx*2.0800},{\dy*1.0365}) + -- ({\dx*2.1000},{\dy*1.0465}) + -- ({\dx*2.1200},{\dy*1.0568}) + -- ({\dx*2.1400},{\dy*1.0675}) + -- ({\dx*2.1600},{\dy*1.0786}) + -- ({\dx*2.1800},{\dy*1.0900}) + -- ({\dx*2.2000},{\dy*1.1018}) + -- ({\dx*2.2200},{\dy*1.1140}) + -- ({\dx*2.2400},{\dy*1.1266}) + -- ({\dx*2.2600},{\dy*1.1395}) + -- ({\dx*2.2800},{\dy*1.1529}) + -- ({\dx*2.3000},{\dy*1.1667}) + -- ({\dx*2.3200},{\dy*1.1809}) + -- ({\dx*2.3400},{\dy*1.1956}) + -- ({\dx*2.3600},{\dy*1.2107}) + -- ({\dx*2.3800},{\dy*1.2262}) + -- ({\dx*2.4000},{\dy*1.2422}) + -- ({\dx*2.4200},{\dy*1.2586}) + -- ({\dx*2.4400},{\dy*1.2756}) + -- ({\dx*2.4600},{\dy*1.2930}) + -- ({\dx*2.4800},{\dy*1.3109}) + -- ({\dx*2.5000},{\dy*1.3293}) + -- ({\dx*2.5200},{\dy*1.3483}) + -- ({\dx*2.5400},{\dy*1.3678}) + -- ({\dx*2.5600},{\dy*1.3878}) + -- ({\dx*2.5800},{\dy*1.4084}) + -- ({\dx*2.6000},{\dy*1.4296}) + -- ({\dx*2.6200},{\dy*1.4514}) + -- ({\dx*2.6400},{\dy*1.4738}) + -- ({\dx*2.6600},{\dy*1.4968}) + -- ({\dx*2.6800},{\dy*1.5204}) + -- ({\dx*2.7000},{\dy*1.5447}) + -- ({\dx*2.7200},{\dy*1.5696}) + -- ({\dx*2.7400},{\dy*1.5953}) + -- ({\dx*2.7600},{\dy*1.6216}) + -- ({\dx*2.7800},{\dy*1.6487}) + -- ({\dx*2.8000},{\dy*1.6765}) + -- ({\dx*2.8200},{\dy*1.7051}) + -- ({\dx*2.8400},{\dy*1.7344}) + -- ({\dx*2.8600},{\dy*1.7646}) + -- ({\dx*2.8800},{\dy*1.7955}) + -- ({\dx*2.9000},{\dy*1.8274}) + -- ({\dx*2.9200},{\dy*1.8600}) + -- ({\dx*2.9400},{\dy*1.8936}) + -- ({\dx*2.9600},{\dy*1.9281}) + -- ({\dx*2.9800},{\dy*1.9636}) + -- ({\dx*3.0000},{\dy*2.0000}) + -- ({\dx*3.0200},{\dy*2.0374}) + -- ({\dx*3.0400},{\dy*2.0759}) + -- ({\dx*3.0600},{\dy*2.1153}) + -- ({\dx*3.0800},{\dy*2.1559}) + -- ({\dx*3.1000},{\dy*2.1976}) + -- ({\dx*3.1200},{\dy*2.2405}) + -- ({\dx*3.1400},{\dy*2.2845}) + -- ({\dx*3.1600},{\dy*2.3297}) + -- ({\dx*3.1800},{\dy*2.3762}) + -- ({\dx*3.2000},{\dy*2.4240}) + -- ({\dx*3.2200},{\dy*2.4731}) + -- ({\dx*3.2400},{\dy*2.5235}) + -- ({\dx*3.2600},{\dy*2.5754}) + -- ({\dx*3.2800},{\dy*2.6287}) + -- ({\dx*3.3000},{\dy*2.6834}) + -- ({\dx*3.3200},{\dy*2.7397}) + -- ({\dx*3.3400},{\dy*2.7976}) + -- ({\dx*3.3600},{\dy*2.8571}) + -- ({\dx*3.3800},{\dy*2.9183}) + -- ({\dx*3.4000},{\dy*2.9812}) + -- ({\dx*3.4200},{\dy*3.0459}) + -- ({\dx*3.4400},{\dy*3.1124}) + -- ({\dx*3.4600},{\dy*3.1807}) + -- ({\dx*3.4800},{\dy*3.2510}) + -- ({\dx*3.5000},{\dy*3.3234}) + -- ({\dx*3.5200},{\dy*3.3977}) + -- ({\dx*3.5400},{\dy*3.4742}) + -- ({\dx*3.5600},{\dy*3.5529}) + -- ({\dx*3.5800},{\dy*3.6338}) + -- ({\dx*3.6000},{\dy*3.7170}) + -- ({\dx*3.6200},{\dy*3.8027}) + -- ({\dx*3.6400},{\dy*3.8908}) + -- ({\dx*3.6600},{\dy*3.9814}) + -- ({\dx*3.6800},{\dy*4.0747}) + -- ({\dx*3.7000},{\dy*4.1707}) + -- ({\dx*3.7200},{\dy*4.2694}) + -- ({\dx*3.7400},{\dy*4.3711}) + -- ({\dx*3.7600},{\dy*4.4757}) + -- ({\dx*3.7800},{\dy*4.5833}) + -- ({\dx*3.8000},{\dy*4.6942}) + -- ({\dx*3.8200},{\dy*4.8083}) + -- ({\dx*3.8400},{\dy*4.9257}) + -- ({\dx*3.8600},{\dy*5.0466}) + -- ({\dx*3.8800},{\dy*5.1711}) + -- ({\dx*3.9000},{\dy*5.2993}) + -- ({\dx*3.9200},{\dy*5.4313}) + -- ({\dx*3.9400},{\dy*5.5673}) + -- ({\dx*3.9600},{\dy*5.7073}) + -- ({\dx*3.9800},{\dy*5.8515}) + -- ({\dx*4.0000},{\dy*6.0000}) + -- ({\dx*4.0200},{\dy*6.1530}) + -- ({\dx*4.0400},{\dy*6.3106}) + -- ({\dx*4.0600},{\dy*6.4730}) + -- ({\dx*4.0800},{\dy*6.6403}) + -- ({\dx*4.0810},{\dy*6.6488})} +\def\gammaone{({\dx*-0.9810},{\dy*-53.0814}) + -- ({\dx*-0.9800},{\dy*-50.4512}) + -- ({\dx*-0.9600},{\dy*-25.4802}) + -- ({\dx*-0.9400},{\dy*-17.1763}) + -- ({\dx*-0.9200},{\dy*-13.0397}) + -- ({\dx*-0.9000},{\dy*-10.5706}) + -- ({\dx*-0.8800},{\dy*-8.9355}) + -- ({\dx*-0.8600},{\dy*-7.7775}) + -- ({\dx*-0.8400},{\dy*-6.9182}) + -- ({\dx*-0.8200},{\dy*-6.2583}) + -- ({\dx*-0.8000},{\dy*-5.7386}) + -- ({\dx*-0.7800},{\dy*-5.3211}) + -- ({\dx*-0.7600},{\dy*-4.9809}) + -- ({\dx*-0.7400},{\dy*-4.7006}) + -- ({\dx*-0.7200},{\dy*-4.4678}) + -- ({\dx*-0.7000},{\dy*-4.2737}) + -- ({\dx*-0.6800},{\dy*-4.1114}) + -- ({\dx*-0.6600},{\dy*-3.9760}) + -- ({\dx*-0.6400},{\dy*-3.8636}) + -- ({\dx*-0.6200},{\dy*-3.7714}) + -- ({\dx*-0.6000},{\dy*-3.6969}) + -- ({\dx*-0.5800},{\dy*-3.6386}) + -- ({\dx*-0.5600},{\dy*-3.5950}) + -- ({\dx*-0.5400},{\dy*-3.5652}) + -- ({\dx*-0.5200},{\dy*-3.5487}) + -- ({\dx*-0.5000},{\dy*-3.5449}) + -- ({\dx*-0.4800},{\dy*-3.5538}) + -- ({\dx*-0.4600},{\dy*-3.5756}) + -- ({\dx*-0.4400},{\dy*-3.6105}) + -- ({\dx*-0.4200},{\dy*-3.6594}) + -- ({\dx*-0.4000},{\dy*-3.7230}) + -- ({\dx*-0.3800},{\dy*-3.8027}) + -- ({\dx*-0.3600},{\dy*-3.9004}) + -- ({\dx*-0.3400},{\dy*-4.0181}) + -- ({\dx*-0.3200},{\dy*-4.1590}) + -- ({\dx*-0.3000},{\dy*-4.3269}) + -- ({\dx*-0.2800},{\dy*-4.5267}) + -- ({\dx*-0.2600},{\dy*-4.7652}) + -- ({\dx*-0.2400},{\dy*-5.0514}) + -- ({\dx*-0.2200},{\dy*-5.3976}) + -- ({\dx*-0.2000},{\dy*-5.8211}) + -- ({\dx*-0.1800},{\dy*-6.3472}) + -- ({\dx*-0.1600},{\dy*-7.0135}) + -- ({\dx*-0.1400},{\dy*-7.8795}) + -- ({\dx*-0.1200},{\dy*-9.0442}) + -- ({\dx*-0.1000},{\dy*-10.6863}) + -- ({\dx*-0.0800},{\dy*-13.1627}) + -- ({\dx*-0.0600},{\dy*-17.3067}) + -- ({\dx*-0.0400},{\dy*-25.6183}) + -- ({\dx*-0.0200},{\dy*-50.5974}) + -- ({\dx*-0.0190},{\dy*-53.2279})} +\def\gammatwo{({\dx*-1.9810},{\dy*26.7952}) + -- ({\dx*-1.9800},{\dy*25.4804}) + -- ({\dx*-1.9600},{\dy*13.0001}) + -- ({\dx*-1.9400},{\dy*8.8538}) + -- ({\dx*-1.9200},{\dy*6.7915}) + -- ({\dx*-1.9000},{\dy*5.5635}) + -- ({\dx*-1.8800},{\dy*4.7529}) + -- ({\dx*-1.8600},{\dy*4.1815}) + -- ({\dx*-1.8400},{\dy*3.7599}) + -- ({\dx*-1.8200},{\dy*3.4386}) + -- ({\dx*-1.8000},{\dy*3.1881}) + -- ({\dx*-1.7800},{\dy*2.9894}) + -- ({\dx*-1.7600},{\dy*2.8301}) + -- ({\dx*-1.7400},{\dy*2.7015}) + -- ({\dx*-1.7200},{\dy*2.5976}) + -- ({\dx*-1.7000},{\dy*2.5139}) + -- ({\dx*-1.6800},{\dy*2.4473}) + -- ({\dx*-1.6600},{\dy*2.3952}) + -- ({\dx*-1.6400},{\dy*2.3559}) + -- ({\dx*-1.6200},{\dy*2.3280}) + -- ({\dx*-1.6000},{\dy*2.3106}) + -- ({\dx*-1.5800},{\dy*2.3029}) + -- ({\dx*-1.5600},{\dy*2.3045}) + -- ({\dx*-1.5400},{\dy*2.3151}) + -- ({\dx*-1.5200},{\dy*2.3346}) + -- ({\dx*-1.5000},{\dy*2.3633}) + -- ({\dx*-1.4800},{\dy*2.4012}) + -- ({\dx*-1.4600},{\dy*2.4490}) + -- ({\dx*-1.4400},{\dy*2.5073}) + -- ({\dx*-1.4200},{\dy*2.5770}) + -- ({\dx*-1.4000},{\dy*2.6593}) + -- ({\dx*-1.3800},{\dy*2.7556}) + -- ({\dx*-1.3600},{\dy*2.8679}) + -- ({\dx*-1.3400},{\dy*2.9986}) + -- ({\dx*-1.3200},{\dy*3.1508}) + -- ({\dx*-1.3000},{\dy*3.3283}) + -- ({\dx*-1.2800},{\dy*3.5365}) + -- ({\dx*-1.2600},{\dy*3.7819}) + -- ({\dx*-1.2400},{\dy*4.0737}) + -- ({\dx*-1.2200},{\dy*4.4243}) + -- ({\dx*-1.2000},{\dy*4.8510}) + -- ({\dx*-1.1800},{\dy*5.3790}) + -- ({\dx*-1.1600},{\dy*6.0461}) + -- ({\dx*-1.1400},{\dy*6.9118}) + -- ({\dx*-1.1200},{\dy*8.0752}) + -- ({\dx*-1.1000},{\dy*9.7148}) + -- ({\dx*-1.0800},{\dy*12.1877}) + -- ({\dx*-1.0600},{\dy*16.3271}) + -- ({\dx*-1.0400},{\dy*24.6330}) + -- ({\dx*-1.0200},{\dy*49.6053}) + -- ({\dx*-1.0190},{\dy*52.2354})} +\def\gammathree{({\dx*-2.9810},{\dy*-8.9887}) + -- ({\dx*-2.9800},{\dy*-8.5505}) + -- ({\dx*-2.9600},{\dy*-4.3919}) + -- ({\dx*-2.9400},{\dy*-3.0115}) + -- ({\dx*-2.9200},{\dy*-2.3259}) + -- ({\dx*-2.9000},{\dy*-1.9184}) + -- ({\dx*-2.8800},{\dy*-1.6503}) + -- ({\dx*-2.8600},{\dy*-1.4621}) + -- ({\dx*-2.8400},{\dy*-1.3239}) + -- ({\dx*-2.8200},{\dy*-1.2194}) + -- ({\dx*-2.8000},{\dy*-1.1386}) + -- ({\dx*-2.7800},{\dy*-1.0753}) + -- ({\dx*-2.7600},{\dy*-1.0254}) + -- ({\dx*-2.7400},{\dy*-0.9859}) + -- ({\dx*-2.7200},{\dy*-0.9550}) + -- ({\dx*-2.7000},{\dy*-0.9311}) + -- ({\dx*-2.6800},{\dy*-0.9132}) + -- ({\dx*-2.6600},{\dy*-0.9004}) + -- ({\dx*-2.6400},{\dy*-0.8924}) + -- ({\dx*-2.6200},{\dy*-0.8886}) + -- ({\dx*-2.6000},{\dy*-0.8887}) + -- ({\dx*-2.5800},{\dy*-0.8926}) + -- ({\dx*-2.5600},{\dy*-0.9002}) + -- ({\dx*-2.5400},{\dy*-0.9115}) + -- ({\dx*-2.5200},{\dy*-0.9264}) + -- ({\dx*-2.5000},{\dy*-0.9453}) + -- ({\dx*-2.4800},{\dy*-0.9682}) + -- ({\dx*-2.4600},{\dy*-0.9955}) + -- ({\dx*-2.4400},{\dy*-1.0276}) + -- ({\dx*-2.4200},{\dy*-1.0649}) + -- ({\dx*-2.4000},{\dy*-1.1080}) + -- ({\dx*-2.3800},{\dy*-1.1578}) + -- ({\dx*-2.3600},{\dy*-1.2152}) + -- ({\dx*-2.3400},{\dy*-1.2815}) + -- ({\dx*-2.3200},{\dy*-1.3581}) + -- ({\dx*-2.3000},{\dy*-1.4471}) + -- ({\dx*-2.2800},{\dy*-1.5511}) + -- ({\dx*-2.2600},{\dy*-1.6734}) + -- ({\dx*-2.2400},{\dy*-1.8186}) + -- ({\dx*-2.2200},{\dy*-1.9929}) + -- ({\dx*-2.2000},{\dy*-2.2050}) + -- ({\dx*-2.1800},{\dy*-2.4674}) + -- ({\dx*-2.1600},{\dy*-2.7991}) + -- ({\dx*-2.1400},{\dy*-3.2298}) + -- ({\dx*-2.1200},{\dy*-3.8091}) + -- ({\dx*-2.1000},{\dy*-4.6261}) + -- ({\dx*-2.0800},{\dy*-5.8595}) + -- ({\dx*-2.0600},{\dy*-7.9258}) + -- ({\dx*-2.0400},{\dy*-12.0750}) + -- ({\dx*-2.0200},{\dy*-24.5571}) + -- ({\dx*-2.0190},{\dy*-25.8719})} +\def\gammafour{({\dx*-3.9950},{\dy*8.3966}) + -- ({\dx*-3.9800},{\dy*2.1484}) + -- ({\dx*-3.9600},{\dy*1.1091}) + -- ({\dx*-3.9400},{\dy*0.7643}) + -- ({\dx*-3.9200},{\dy*0.5933}) + -- ({\dx*-3.9000},{\dy*0.4919}) + -- ({\dx*-3.8800},{\dy*0.4253}) + -- ({\dx*-3.8600},{\dy*0.3788}) + -- ({\dx*-3.8400},{\dy*0.3448}) + -- ({\dx*-3.8200},{\dy*0.3192}) + -- ({\dx*-3.8000},{\dy*0.2996}) + -- ({\dx*-3.7800},{\dy*0.2845}) + -- ({\dx*-3.7600},{\dy*0.2727}) + -- ({\dx*-3.7400},{\dy*0.2636}) + -- ({\dx*-3.7200},{\dy*0.2567}) + -- ({\dx*-3.7000},{\dy*0.2516}) + -- ({\dx*-3.6800},{\dy*0.2481}) + -- ({\dx*-3.6600},{\dy*0.2460}) + -- ({\dx*-3.6400},{\dy*0.2452}) + -- ({\dx*-3.6200},{\dy*0.2455}) + -- ({\dx*-3.6000},{\dy*0.2469}) + -- ({\dx*-3.5800},{\dy*0.2493}) + -- ({\dx*-3.5600},{\dy*0.2529}) + -- ({\dx*-3.5400},{\dy*0.2575}) + -- ({\dx*-3.5200},{\dy*0.2632}) + -- ({\dx*-3.5000},{\dy*0.2701}) + -- ({\dx*-3.4800},{\dy*0.2782}) + -- ({\dx*-3.4600},{\dy*0.2877}) + -- ({\dx*-3.4400},{\dy*0.2987}) + -- ({\dx*-3.4200},{\dy*0.3114}) + -- ({\dx*-3.4000},{\dy*0.3259}) + -- ({\dx*-3.3800},{\dy*0.3425}) + -- ({\dx*-3.3600},{\dy*0.3617}) + -- ({\dx*-3.3400},{\dy*0.3837}) + -- ({\dx*-3.3200},{\dy*0.4091}) + -- ({\dx*-3.3000},{\dy*0.4385}) + -- ({\dx*-3.2800},{\dy*0.4729}) + -- ({\dx*-3.2600},{\dy*0.5133}) + -- ({\dx*-3.2400},{\dy*0.5613}) + -- ({\dx*-3.2200},{\dy*0.6189}) + -- ({\dx*-3.2000},{\dy*0.6891}) + -- ({\dx*-3.1800},{\dy*0.7759}) + -- ({\dx*-3.1600},{\dy*0.8858}) + -- ({\dx*-3.1400},{\dy*1.0286}) + -- ({\dx*-3.1200},{\dy*1.2209}) + -- ({\dx*-3.1000},{\dy*1.4923}) + -- ({\dx*-3.0800},{\dy*1.9024}) + -- ({\dx*-3.0600},{\dy*2.5901}) + -- ({\dx*-3.0400},{\dy*3.9720}) + -- ({\dx*-3.0200},{\dy*8.1315}) + -- ({\dx*-3.0050},{\dy*33.1259})} +\def\gammafive{({\dx*-4.9990},{\dy*-8.3476}) + -- ({\dx*-4.9800},{\dy*-0.4314}) + -- ({\dx*-4.9600},{\dy*-0.2236}) + -- ({\dx*-4.9400},{\dy*-0.1547}) + -- ({\dx*-4.9200},{\dy*-0.1206}) + -- ({\dx*-4.9000},{\dy*-0.1004}) + -- ({\dx*-4.8800},{\dy*-0.0872}) + -- ({\dx*-4.8600},{\dy*-0.0779}) + -- ({\dx*-4.8400},{\dy*-0.0712}) + -- ({\dx*-4.8200},{\dy*-0.0662}) + -- ({\dx*-4.8000},{\dy*-0.0624}) + -- ({\dx*-4.7800},{\dy*-0.0595}) + -- ({\dx*-4.7600},{\dy*-0.0573}) + -- ({\dx*-4.7400},{\dy*-0.0556}) + -- ({\dx*-4.7200},{\dy*-0.0544}) + -- ({\dx*-4.7000},{\dy*-0.0535}) + -- ({\dx*-4.6800},{\dy*-0.0530}) + -- ({\dx*-4.6600},{\dy*-0.0528}) + -- ({\dx*-4.6400},{\dy*-0.0528}) + -- ({\dx*-4.6200},{\dy*-0.0531}) + -- ({\dx*-4.6000},{\dy*-0.0537}) + -- ({\dx*-4.5800},{\dy*-0.0544}) + -- ({\dx*-4.5600},{\dy*-0.0555}) + -- ({\dx*-4.5400},{\dy*-0.0567}) + -- ({\dx*-4.5200},{\dy*-0.0582}) + -- ({\dx*-4.5000},{\dy*-0.0600}) + -- ({\dx*-4.4800},{\dy*-0.0621}) + -- ({\dx*-4.4600},{\dy*-0.0645}) + -- ({\dx*-4.4400},{\dy*-0.0673}) + -- ({\dx*-4.4200},{\dy*-0.0704}) + -- ({\dx*-4.4000},{\dy*-0.0741}) + -- ({\dx*-4.3800},{\dy*-0.0782}) + -- ({\dx*-4.3600},{\dy*-0.0830}) + -- ({\dx*-4.3400},{\dy*-0.0884}) + -- ({\dx*-4.3200},{\dy*-0.0947}) + -- ({\dx*-4.3000},{\dy*-0.1020}) + -- ({\dx*-4.2800},{\dy*-0.1105}) + -- ({\dx*-4.2600},{\dy*-0.1205}) + -- ({\dx*-4.2400},{\dy*-0.1324}) + -- ({\dx*-4.2200},{\dy*-0.1467}) + -- ({\dx*-4.2000},{\dy*-0.1641}) + -- ({\dx*-4.1800},{\dy*-0.1856}) + -- ({\dx*-4.1600},{\dy*-0.2129}) + -- ({\dx*-4.1400},{\dy*-0.2485}) + -- ({\dx*-4.1200},{\dy*-0.2963}) + -- ({\dx*-4.1000},{\dy*-0.3640}) + -- ({\dx*-4.0800},{\dy*-0.4663}) + -- ({\dx*-4.0600},{\dy*-0.6380}) + -- ({\dx*-4.0400},{\dy*-0.9832}) + -- ({\dx*-4.0200},{\dy*-2.0228}) + -- ({\dx*-4.0010},{\dy*-41.6040})} +\def\gammasix{({\dx*-5.9998},{\dy*6.9470}) + -- ({\dx*-5.9800},{\dy*0.0721}) + -- ({\dx*-5.9600},{\dy*0.0375}) + -- ({\dx*-5.9400},{\dy*0.0260}) + -- ({\dx*-5.9200},{\dy*0.0204}) + -- ({\dx*-5.9000},{\dy*0.0170}) + -- ({\dx*-5.8800},{\dy*0.0148}) + -- ({\dx*-5.8600},{\dy*0.0133}) + -- ({\dx*-5.8400},{\dy*0.0122}) + -- ({\dx*-5.8200},{\dy*0.0114}) + -- ({\dx*-5.8000},{\dy*0.0108}) + -- ({\dx*-5.7800},{\dy*0.0103}) + -- ({\dx*-5.7600},{\dy*0.0099}) + -- ({\dx*-5.7400},{\dy*0.0097}) + -- ({\dx*-5.7200},{\dy*0.0095}) + -- ({\dx*-5.7000},{\dy*0.0094}) + -- ({\dx*-5.6800},{\dy*0.0093}) + -- ({\dx*-5.6600},{\dy*0.0093}) + -- ({\dx*-5.6400},{\dy*0.0094}) + -- ({\dx*-5.6200},{\dy*0.0095}) + -- ({\dx*-5.6000},{\dy*0.0096}) + -- ({\dx*-5.5800},{\dy*0.0098}) + -- ({\dx*-5.5600},{\dy*0.0100}) + -- ({\dx*-5.5400},{\dy*0.0102}) + -- ({\dx*-5.5200},{\dy*0.0105}) + -- ({\dx*-5.5000},{\dy*0.0109}) + -- ({\dx*-5.4800},{\dy*0.0113}) + -- ({\dx*-5.4600},{\dy*0.0118}) + -- ({\dx*-5.4400},{\dy*0.0124}) + -- ({\dx*-5.4200},{\dy*0.0130}) + -- ({\dx*-5.4000},{\dy*0.0137}) + -- ({\dx*-5.3800},{\dy*0.0145}) + -- ({\dx*-5.3600},{\dy*0.0155}) + -- ({\dx*-5.3400},{\dy*0.0166}) + -- ({\dx*-5.3200},{\dy*0.0178}) + -- ({\dx*-5.3000},{\dy*0.0192}) + -- ({\dx*-5.2800},{\dy*0.0209}) + -- ({\dx*-5.2600},{\dy*0.0229}) + -- ({\dx*-5.2400},{\dy*0.0253}) + -- ({\dx*-5.2200},{\dy*0.0281}) + -- ({\dx*-5.2000},{\dy*0.0316}) + -- ({\dx*-5.1800},{\dy*0.0358}) + -- ({\dx*-5.1600},{\dy*0.0413}) + -- ({\dx*-5.1400},{\dy*0.0483}) + -- ({\dx*-5.1200},{\dy*0.0579}) + -- ({\dx*-5.1000},{\dy*0.0714}) + -- ({\dx*-5.0800},{\dy*0.0918}) + -- ({\dx*-5.0600},{\dy*0.1261}) + -- ({\dx*-5.0400},{\dy*0.1951}) + -- ({\dx*-5.0200},{\dy*0.4029}) + -- ({\dx*-5.0002},{\dy*41.6525})} +\def\gammasinplus{({\dx*0.0190},{\dy*52.1325}) + -- ({\dx*0.0200},{\dy*49.5050}) + -- ({\dx*0.0400},{\dy*24.5863}) + -- ({\dx*0.0600},{\dy*16.3331}) + -- ({\dx*0.0800},{\dy*12.2453}) + -- ({\dx*0.1000},{\dy*9.8225}) + -- ({\dx*0.1200},{\dy*8.2314}) + -- ({\dx*0.1400},{\dy*7.1145}) + -- ({\dx*0.1600},{\dy*6.2930}) + -- ({\dx*0.1800},{\dy*5.6676}) + -- ({\dx*0.2000},{\dy*5.1786}) + -- ({\dx*0.2200},{\dy*4.7879}) + -- ({\dx*0.2400},{\dy*4.4701}) + -- ({\dx*0.2600},{\dy*4.2074}) + -- ({\dx*0.2800},{\dy*3.9874}) + -- ({\dx*0.3000},{\dy*3.8006}) + -- ({\dx*0.3200},{\dy*3.6401}) + -- ({\dx*0.3400},{\dy*3.5005}) + -- ({\dx*0.3600},{\dy*3.3776}) + -- ({\dx*0.3800},{\dy*3.2680}) + -- ({\dx*0.4000},{\dy*3.1692}) + -- ({\dx*0.4200},{\dy*3.0790}) + -- ({\dx*0.4400},{\dy*2.9955}) + -- ({\dx*0.4600},{\dy*2.9173}) + -- ({\dx*0.4800},{\dy*2.8433}) + -- ({\dx*0.5000},{\dy*2.7725}) + -- ({\dx*0.5200},{\dy*2.7039}) + -- ({\dx*0.5400},{\dy*2.6369}) + -- ({\dx*0.5600},{\dy*2.5709}) + -- ({\dx*0.5800},{\dy*2.5055}) + -- ({\dx*0.6000},{\dy*2.4402}) + -- ({\dx*0.6200},{\dy*2.3748}) + -- ({\dx*0.6400},{\dy*2.3090}) + -- ({\dx*0.6600},{\dy*2.2425}) + -- ({\dx*0.6800},{\dy*2.1752}) + -- ({\dx*0.7000},{\dy*2.1071}) + -- ({\dx*0.7200},{\dy*2.0380}) + -- ({\dx*0.7400},{\dy*1.9679}) + -- ({\dx*0.7600},{\dy*1.8969}) + -- ({\dx*0.7800},{\dy*1.8249}) + -- ({\dx*0.8000},{\dy*1.7520}) + -- ({\dx*0.8200},{\dy*1.6783}) + -- ({\dx*0.8400},{\dy*1.6039}) + -- ({\dx*0.8600},{\dy*1.5289}) + -- ({\dx*0.8800},{\dy*1.4534}) + -- ({\dx*0.9000},{\dy*1.3776}) + -- ({\dx*0.9200},{\dy*1.3017}) + -- ({\dx*0.9400},{\dy*1.2258}) + -- ({\dx*0.9600},{\dy*1.1501}) + -- ({\dx*0.9800},{\dy*1.0747}) + -- ({\dx*1.0000},{\dy*1.0000}) + -- ({\dx*1.0200},{\dy*0.9261}) + -- ({\dx*1.0400},{\dy*0.8531}) + -- ({\dx*1.0600},{\dy*0.7814}) + -- ({\dx*1.0800},{\dy*0.7110}) + -- ({\dx*1.1000},{\dy*0.6423}) + -- ({\dx*1.1200},{\dy*0.5755}) + -- ({\dx*1.1400},{\dy*0.5106}) + -- ({\dx*1.1600},{\dy*0.4480}) + -- ({\dx*1.1800},{\dy*0.3879}) + -- ({\dx*1.2000},{\dy*0.3304}) + -- ({\dx*1.2200},{\dy*0.2757}) + -- ({\dx*1.2400},{\dy*0.2240}) + -- ({\dx*1.2600},{\dy*0.1754}) + -- ({\dx*1.2800},{\dy*0.1302}) + -- ({\dx*1.3000},{\dy*0.0885}) + -- ({\dx*1.3200},{\dy*0.0503}) + -- ({\dx*1.3400},{\dy*0.0159}) + -- ({\dx*1.3600},{\dy*-0.0146}) + -- ({\dx*1.3800},{\dy*-0.0412}) + -- ({\dx*1.4000},{\dy*-0.0638}) + -- ({\dx*1.4200},{\dy*-0.0822}) + -- ({\dx*1.4400},{\dy*-0.0965}) + -- ({\dx*1.4600},{\dy*-0.1065}) + -- ({\dx*1.4800},{\dy*-0.1123}) + -- ({\dx*1.5000},{\dy*-0.1138}) + -- ({\dx*1.5200},{\dy*-0.1110}) + -- ({\dx*1.5400},{\dy*-0.1039}) + -- ({\dx*1.5600},{\dy*-0.0926}) + -- ({\dx*1.5800},{\dy*-0.0772}) + -- ({\dx*1.6000},{\dy*-0.0575}) + -- ({\dx*1.6200},{\dy*-0.0339}) + -- ({\dx*1.6400},{\dy*-0.0062}) + -- ({\dx*1.6600},{\dy*0.0254}) + -- ({\dx*1.6800},{\dy*0.0607}) + -- ({\dx*1.7000},{\dy*0.0996}) + -- ({\dx*1.7200},{\dy*0.1421}) + -- ({\dx*1.7400},{\dy*0.1879}) + -- ({\dx*1.7600},{\dy*0.2368}) + -- ({\dx*1.7800},{\dy*0.2888}) + -- ({\dx*1.8000},{\dy*0.3436}) + -- ({\dx*1.8200},{\dy*0.4010}) + -- ({\dx*1.8400},{\dy*0.4609}) + -- ({\dx*1.8600},{\dy*0.5229}) + -- ({\dx*1.8800},{\dy*0.5869}) + -- ({\dx*1.9000},{\dy*0.6527}) + -- ({\dx*1.9200},{\dy*0.7201}) + -- ({\dx*1.9400},{\dy*0.7887}) + -- ({\dx*1.9600},{\dy*0.8584}) + -- ({\dx*1.9800},{\dy*0.9289}) + -- ({\dx*2.0000},{\dy*1.0000}) + -- ({\dx*2.0200},{\dy*1.0714}) + -- ({\dx*2.0400},{\dy*1.1429}) + -- ({\dx*2.0600},{\dy*1.2142}) + -- ({\dx*2.0800},{\dy*1.2852}) + -- ({\dx*2.1000},{\dy*1.3555}) + -- ({\dx*2.1200},{\dy*1.4249}) + -- ({\dx*2.1400},{\dy*1.4933}) + -- ({\dx*2.1600},{\dy*1.5603}) + -- ({\dx*2.1800},{\dy*1.6258}) + -- ({\dx*2.2000},{\dy*1.6896}) + -- ({\dx*2.2200},{\dy*1.7514}) + -- ({\dx*2.2400},{\dy*1.8111}) + -- ({\dx*2.2600},{\dy*1.8685}) + -- ({\dx*2.2800},{\dy*1.9234}) + -- ({\dx*2.3000},{\dy*1.9757}) + -- ({\dx*2.3200},{\dy*2.0253}) + -- ({\dx*2.3400},{\dy*2.0719}) + -- ({\dx*2.3600},{\dy*2.1155}) + -- ({\dx*2.3800},{\dy*2.1560}) + -- ({\dx*2.4000},{\dy*2.1932}) + -- ({\dx*2.4200},{\dy*2.2272}) + -- ({\dx*2.4400},{\dy*2.2578}) + -- ({\dx*2.4600},{\dy*2.2851}) + -- ({\dx*2.4800},{\dy*2.3089}) + -- ({\dx*2.5000},{\dy*2.3293}) + -- ({\dx*2.5200},{\dy*2.3463}) + -- ({\dx*2.5400},{\dy*2.3599}) + -- ({\dx*2.5600},{\dy*2.3701}) + -- ({\dx*2.5800},{\dy*2.3770}) + -- ({\dx*2.6000},{\dy*2.3807}) + -- ({\dx*2.6200},{\dy*2.3812}) + -- ({\dx*2.6400},{\dy*2.3786}) + -- ({\dx*2.6600},{\dy*2.3731}) + -- ({\dx*2.6800},{\dy*2.3647}) + -- ({\dx*2.7000},{\dy*2.3537}) + -- ({\dx*2.7200},{\dy*2.3402}) + -- ({\dx*2.7400},{\dy*2.3242}) + -- ({\dx*2.7600},{\dy*2.3062}) + -- ({\dx*2.7800},{\dy*2.2861}) + -- ({\dx*2.8000},{\dy*2.2643}) + -- ({\dx*2.8200},{\dy*2.2409}) + -- ({\dx*2.8400},{\dy*2.2162}) + -- ({\dx*2.8600},{\dy*2.1903}) + -- ({\dx*2.8800},{\dy*2.1637}) + -- ({\dx*2.9000},{\dy*2.1364}) + -- ({\dx*2.9200},{\dy*2.1087}) + -- ({\dx*2.9400},{\dy*2.0810}) + -- ({\dx*2.9600},{\dy*2.0535}) + -- ({\dx*2.9800},{\dy*2.0264}) + -- ({\dx*3.0000},{\dy*2.0000}) + -- ({\dx*3.0200},{\dy*1.9746}) + -- ({\dx*3.0400},{\dy*1.9505}) + -- ({\dx*3.0600},{\dy*1.9280}) + -- ({\dx*3.0800},{\dy*1.9072}) + -- ({\dx*3.1000},{\dy*1.8886}) + -- ({\dx*3.1200},{\dy*1.8723}) + -- ({\dx*3.1400},{\dy*1.8587}) + -- ({\dx*3.1600},{\dy*1.8480}) + -- ({\dx*3.1800},{\dy*1.8404}) + -- ({\dx*3.2000},{\dy*1.8362}) + -- ({\dx*3.2200},{\dy*1.8356}) + -- ({\dx*3.2400},{\dy*1.8390}) + -- ({\dx*3.2600},{\dy*1.8464}) + -- ({\dx*3.2800},{\dy*1.8581}) + -- ({\dx*3.3000},{\dy*1.8744}) + -- ({\dx*3.3200},{\dy*1.8954}) + -- ({\dx*3.3400},{\dy*1.9213}) + -- ({\dx*3.3600},{\dy*1.9523}) + -- ({\dx*3.3800},{\dy*1.9885}) + -- ({\dx*3.4000},{\dy*2.0301}) + -- ({\dx*3.4200},{\dy*2.0773}) + -- ({\dx*3.4400},{\dy*2.1301}) + -- ({\dx*3.4600},{\dy*2.1886}) + -- ({\dx*3.4800},{\dy*2.2530}) + -- ({\dx*3.5000},{\dy*2.3234}) + -- ({\dx*3.5200},{\dy*2.3997}) + -- ({\dx*3.5400},{\dy*2.4821}) + -- ({\dx*3.5600},{\dy*2.5706}) + -- ({\dx*3.5800},{\dy*2.6652}) + -- ({\dx*3.6000},{\dy*2.7660}) + -- ({\dx*3.6200},{\dy*2.8729}) + -- ({\dx*3.6400},{\dy*2.9859}) + -- ({\dx*3.6600},{\dy*3.1051}) + -- ({\dx*3.6800},{\dy*3.2303}) + -- ({\dx*3.7000},{\dy*3.3616}) + -- ({\dx*3.7200},{\dy*3.4989}) + -- ({\dx*3.7400},{\dy*3.6421}) + -- ({\dx*3.7600},{\dy*3.7911}) + -- ({\dx*3.7800},{\dy*3.9459}) + -- ({\dx*3.8000},{\dy*4.1064}) + -- ({\dx*3.8200},{\dy*4.2724}) + -- ({\dx*3.8400},{\dy*4.4440}) + -- ({\dx*3.8600},{\dy*4.6209}) + -- ({\dx*3.8800},{\dy*4.8030}) + -- ({\dx*3.9000},{\dy*4.9903}) + -- ({\dx*3.9200},{\dy*5.1826}) + -- ({\dx*3.9400},{\dy*5.3799}) + -- ({\dx*3.9600},{\dy*5.5819}) + -- ({\dx*3.9800},{\dy*5.7887}) + -- ({\dx*4.0000},{\dy*6.0000}) + -- ({\dx*4.0200},{\dy*6.2158}) + -- ({\dx*4.0400},{\dy*6.4359}) + -- ({\dx*4.0600},{\dy*6.6603}) + -- ({\dx*4.0800},{\dy*6.8889}) + -- ({\dx*4.0810},{\dy*6.9005})} +\def\gammasinone{({\dx*-0.9810},{\dy*-53.1410}) + -- ({\dx*-0.9800},{\dy*-50.5140}) + -- ({\dx*-0.9600},{\dy*-25.6055}) + -- ({\dx*-0.9400},{\dy*-17.3637}) + -- ({\dx*-0.9200},{\dy*-13.2884}) + -- ({\dx*-0.9000},{\dy*-10.8796}) + -- ({\dx*-0.8800},{\dy*-9.3036}) + -- ({\dx*-0.8600},{\dy*-8.2033}) + -- ({\dx*-0.8400},{\dy*-7.3999}) + -- ({\dx*-0.8200},{\dy*-6.7941}) + -- ({\dx*-0.8000},{\dy*-6.3263}) + -- ({\dx*-0.7800},{\dy*-5.9586}) + -- ({\dx*-0.7600},{\dy*-5.6655}) + -- ({\dx*-0.7400},{\dy*-5.4296}) + -- ({\dx*-0.7200},{\dy*-5.2384}) + -- ({\dx*-0.7000},{\dy*-5.0827}) + -- ({\dx*-0.6800},{\dy*-4.9557}) + -- ({\dx*-0.6600},{\dy*-4.8523}) + -- ({\dx*-0.6400},{\dy*-4.7685}) + -- ({\dx*-0.6200},{\dy*-4.7012}) + -- ({\dx*-0.6000},{\dy*-4.6480}) + -- ({\dx*-0.5800},{\dy*-4.6072}) + -- ({\dx*-0.5600},{\dy*-4.5773}) + -- ({\dx*-0.5400},{\dy*-4.5573}) + -- ({\dx*-0.5200},{\dy*-4.5467}) + -- ({\dx*-0.5000},{\dy*-4.5449}) + -- ({\dx*-0.4800},{\dy*-4.5519}) + -- ({\dx*-0.4600},{\dy*-4.5677}) + -- ({\dx*-0.4400},{\dy*-4.5928}) + -- ({\dx*-0.4200},{\dy*-4.6279}) + -- ({\dx*-0.4000},{\dy*-4.6740}) + -- ({\dx*-0.3800},{\dy*-4.7325}) + -- ({\dx*-0.3600},{\dy*-4.8052}) + -- ({\dx*-0.3400},{\dy*-4.8944}) + -- ({\dx*-0.3200},{\dy*-5.0033}) + -- ({\dx*-0.3000},{\dy*-5.1359}) + -- ({\dx*-0.2800},{\dy*-5.2972}) + -- ({\dx*-0.2600},{\dy*-5.4942}) + -- ({\dx*-0.2400},{\dy*-5.7359}) + -- ({\dx*-0.2200},{\dy*-6.0350}) + -- ({\dx*-0.2000},{\dy*-6.4089}) + -- ({\dx*-0.1800},{\dy*-6.8830}) + -- ({\dx*-0.1600},{\dy*-7.4952}) + -- ({\dx*-0.1400},{\dy*-8.3052}) + -- ({\dx*-0.1200},{\dy*-9.4124}) + -- ({\dx*-0.1000},{\dy*-10.9953}) + -- ({\dx*-0.0800},{\dy*-13.4114}) + -- ({\dx*-0.0600},{\dy*-17.4941}) + -- ({\dx*-0.0400},{\dy*-25.7436}) + -- ({\dx*-0.0200},{\dy*-50.6602}) + -- ({\dx*-0.0190},{\dy*-53.2876})} +\def\gammasintwo{({\dx*-1.9810},{\dy*26.8549}) + -- ({\dx*-1.9800},{\dy*25.5432}) + -- ({\dx*-1.9600},{\dy*13.1254}) + -- ({\dx*-1.9400},{\dy*9.0411}) + -- ({\dx*-1.9200},{\dy*7.0402}) + -- ({\dx*-1.9000},{\dy*5.8725}) + -- ({\dx*-1.8800},{\dy*5.1211}) + -- ({\dx*-1.8600},{\dy*4.6073}) + -- ({\dx*-1.8400},{\dy*4.2416}) + -- ({\dx*-1.8200},{\dy*3.9745}) + -- ({\dx*-1.8000},{\dy*3.7759}) + -- ({\dx*-1.7800},{\dy*3.6268}) + -- ({\dx*-1.7600},{\dy*3.5146}) + -- ({\dx*-1.7400},{\dy*3.4305}) + -- ({\dx*-1.7200},{\dy*3.3681}) + -- ({\dx*-1.7000},{\dy*3.3229}) + -- ({\dx*-1.6800},{\dy*3.2916}) + -- ({\dx*-1.6600},{\dy*3.2715}) + -- ({\dx*-1.6400},{\dy*3.2607}) + -- ({\dx*-1.6200},{\dy*3.2578}) + -- ({\dx*-1.6000},{\dy*3.2616}) + -- ({\dx*-1.5800},{\dy*3.2715}) + -- ({\dx*-1.5600},{\dy*3.2868}) + -- ({\dx*-1.5400},{\dy*3.3072}) + -- ({\dx*-1.5200},{\dy*3.3327}) + -- ({\dx*-1.5000},{\dy*3.3633}) + -- ({\dx*-1.4800},{\dy*3.3993}) + -- ({\dx*-1.4600},{\dy*3.4412}) + -- ({\dx*-1.4400},{\dy*3.4896}) + -- ({\dx*-1.4200},{\dy*3.5456}) + -- ({\dx*-1.4000},{\dy*3.6103}) + -- ({\dx*-1.3800},{\dy*3.6854}) + -- ({\dx*-1.3600},{\dy*3.7727}) + -- ({\dx*-1.3400},{\dy*3.8749}) + -- ({\dx*-1.3200},{\dy*3.9951}) + -- ({\dx*-1.3000},{\dy*4.1374}) + -- ({\dx*-1.2800},{\dy*4.3070}) + -- ({\dx*-1.2600},{\dy*4.5109}) + -- ({\dx*-1.2400},{\dy*4.7583}) + -- ({\dx*-1.2200},{\dy*5.0617}) + -- ({\dx*-1.2000},{\dy*5.4387}) + -- ({\dx*-1.1800},{\dy*5.9148}) + -- ({\dx*-1.1600},{\dy*6.5279}) + -- ({\dx*-1.1400},{\dy*7.3376}) + -- ({\dx*-1.1200},{\dy*8.4433}) + -- ({\dx*-1.1000},{\dy*10.0238}) + -- ({\dx*-1.0800},{\dy*12.4364}) + -- ({\dx*-1.0600},{\dy*16.5145}) + -- ({\dx*-1.0400},{\dy*24.7583}) + -- ({\dx*-1.0200},{\dy*49.6681}) + -- ({\dx*-1.0190},{\dy*52.2951})} +\def\gammasinthree{({\dx*-2.9810},{\dy*-9.0483}) + -- ({\dx*-2.9800},{\dy*-8.6133}) + -- ({\dx*-2.9600},{\dy*-4.5173}) + -- ({\dx*-2.9400},{\dy*-3.1989}) + -- ({\dx*-2.9200},{\dy*-2.5746}) + -- ({\dx*-2.9000},{\dy*-2.2274}) + -- ({\dx*-2.8800},{\dy*-2.0184}) + -- ({\dx*-2.8600},{\dy*-1.8878}) + -- ({\dx*-2.8400},{\dy*-1.8057}) + -- ({\dx*-2.8200},{\dy*-1.7552}) + -- ({\dx*-2.8000},{\dy*-1.7264}) + -- ({\dx*-2.7800},{\dy*-1.7127}) + -- ({\dx*-2.7600},{\dy*-1.7099}) + -- ({\dx*-2.7400},{\dy*-1.7149}) + -- ({\dx*-2.7200},{\dy*-1.7255}) + -- ({\dx*-2.7000},{\dy*-1.7401}) + -- ({\dx*-2.6800},{\dy*-1.7575}) + -- ({\dx*-2.6600},{\dy*-1.7768}) + -- ({\dx*-2.6400},{\dy*-1.7972}) + -- ({\dx*-2.6200},{\dy*-1.8183}) + -- ({\dx*-2.6000},{\dy*-1.8397}) + -- ({\dx*-2.5800},{\dy*-1.8612}) + -- ({\dx*-2.5600},{\dy*-1.8825}) + -- ({\dx*-2.5400},{\dy*-1.9036}) + -- ({\dx*-2.5200},{\dy*-1.9245}) + -- ({\dx*-2.5000},{\dy*-1.9453}) + -- ({\dx*-2.4800},{\dy*-1.9663}) + -- ({\dx*-2.4600},{\dy*-1.9877}) + -- ({\dx*-2.4400},{\dy*-2.0099}) + -- ({\dx*-2.4200},{\dy*-2.0335}) + -- ({\dx*-2.4000},{\dy*-2.0591}) + -- ({\dx*-2.3800},{\dy*-2.0876}) + -- ({\dx*-2.3600},{\dy*-2.1200}) + -- ({\dx*-2.3400},{\dy*-2.1578}) + -- ({\dx*-2.3200},{\dy*-2.2024}) + -- ({\dx*-2.3000},{\dy*-2.2561}) + -- ({\dx*-2.2800},{\dy*-2.3216}) + -- ({\dx*-2.2600},{\dy*-2.4024}) + -- ({\dx*-2.2400},{\dy*-2.5032}) + -- ({\dx*-2.2200},{\dy*-2.6303}) + -- ({\dx*-2.2000},{\dy*-2.7928}) + -- ({\dx*-2.1800},{\dy*-3.0032}) + -- ({\dx*-2.1600},{\dy*-3.2809}) + -- ({\dx*-2.1400},{\dy*-3.6556}) + -- ({\dx*-2.1200},{\dy*-4.1772}) + -- ({\dx*-2.1000},{\dy*-4.9351}) + -- ({\dx*-2.0800},{\dy*-6.1082}) + -- ({\dx*-2.0600},{\dy*-8.1132}) + -- ({\dx*-2.0400},{\dy*-12.2003}) + -- ({\dx*-2.0200},{\dy*-24.6199}) + -- ({\dx*-2.0190},{\dy*-25.9316})} +\def\gammasinfour{({\dx*-3.9950},{\dy*8.4124}) + -- ({\dx*-3.9800},{\dy*2.2112}) + -- ({\dx*-3.9600},{\dy*1.2344}) + -- ({\dx*-3.9400},{\dy*0.9517}) + -- ({\dx*-3.9200},{\dy*0.8420}) + -- ({\dx*-3.9000},{\dy*0.8009}) + -- ({\dx*-3.8800},{\dy*0.7935}) + -- ({\dx*-3.8600},{\dy*0.8045}) + -- ({\dx*-3.8400},{\dy*0.8265}) + -- ({\dx*-3.8200},{\dy*0.8550}) + -- ({\dx*-3.8000},{\dy*0.8874}) + -- ({\dx*-3.7800},{\dy*0.9219}) + -- ({\dx*-3.7600},{\dy*0.9573}) + -- ({\dx*-3.7400},{\dy*0.9926}) + -- ({\dx*-3.7200},{\dy*1.0272}) + -- ({\dx*-3.7000},{\dy*1.0607}) + -- ({\dx*-3.6800},{\dy*1.0925}) + -- ({\dx*-3.6600},{\dy*1.1223}) + -- ({\dx*-3.6400},{\dy*1.1500}) + -- ({\dx*-3.6200},{\dy*1.1752}) + -- ({\dx*-3.6000},{\dy*1.1979}) + -- ({\dx*-3.5800},{\dy*1.2179}) + -- ({\dx*-3.5600},{\dy*1.2351}) + -- ({\dx*-3.5400},{\dy*1.2496}) + -- ({\dx*-3.5200},{\dy*1.2612}) + -- ({\dx*-3.5000},{\dy*1.2701}) + -- ({\dx*-3.4800},{\dy*1.2763}) + -- ({\dx*-3.4600},{\dy*1.2798}) + -- ({\dx*-3.4400},{\dy*1.2810}) + -- ({\dx*-3.4200},{\dy*1.2800}) + -- ({\dx*-3.4000},{\dy*1.2769}) + -- ({\dx*-3.3800},{\dy*1.2723}) + -- ({\dx*-3.3600},{\dy*1.2665}) + -- ({\dx*-3.3400},{\dy*1.2600}) + -- ({\dx*-3.3200},{\dy*1.2534}) + -- ({\dx*-3.3000},{\dy*1.2475}) + -- ({\dx*-3.2800},{\dy*1.2434}) + -- ({\dx*-3.2600},{\dy*1.2423}) + -- ({\dx*-3.2400},{\dy*1.2458}) + -- ({\dx*-3.2200},{\dy*1.2563}) + -- ({\dx*-3.2000},{\dy*1.2768}) + -- ({\dx*-3.1800},{\dy*1.3117}) + -- ({\dx*-3.1600},{\dy*1.3676}) + -- ({\dx*-3.1400},{\dy*1.4544}) + -- ({\dx*-3.1200},{\dy*1.5890}) + -- ({\dx*-3.1000},{\dy*1.8013}) + -- ({\dx*-3.0800},{\dy*2.1511}) + -- ({\dx*-3.0600},{\dy*2.7775}) + -- ({\dx*-3.0400},{\dy*4.0974}) + -- ({\dx*-3.0200},{\dy*8.1943}) + -- ({\dx*-3.0050},{\dy*33.1416})} +\def\gammasinfive{({\dx*-4.9990},{\dy*-8.3507}) + -- ({\dx*-4.9800},{\dy*-0.4942}) + -- ({\dx*-4.9600},{\dy*-0.3489}) + -- ({\dx*-4.9400},{\dy*-0.3421}) + -- ({\dx*-4.9200},{\dy*-0.3693}) + -- ({\dx*-4.9000},{\dy*-0.4094}) + -- ({\dx*-4.8800},{\dy*-0.4553}) + -- ({\dx*-4.8600},{\dy*-0.5037}) + -- ({\dx*-4.8400},{\dy*-0.5530}) + -- ({\dx*-4.8200},{\dy*-0.6021}) + -- ({\dx*-4.8000},{\dy*-0.6502}) + -- ({\dx*-4.7800},{\dy*-0.6969}) + -- ({\dx*-4.7600},{\dy*-0.7418}) + -- ({\dx*-4.7400},{\dy*-0.7846}) + -- ({\dx*-4.7200},{\dy*-0.8249}) + -- ({\dx*-4.7000},{\dy*-0.8626}) + -- ({\dx*-4.6800},{\dy*-0.8973}) + -- ({\dx*-4.6600},{\dy*-0.9291}) + -- ({\dx*-4.6400},{\dy*-0.9577}) + -- ({\dx*-4.6200},{\dy*-0.9829}) + -- ({\dx*-4.6000},{\dy*-1.0047}) + -- ({\dx*-4.5800},{\dy*-1.0230}) + -- ({\dx*-4.5600},{\dy*-1.0377}) + -- ({\dx*-4.5400},{\dy*-1.0488}) + -- ({\dx*-4.5200},{\dy*-1.0563}) + -- ({\dx*-4.5000},{\dy*-1.0600}) + -- ({\dx*-4.4800},{\dy*-1.0601}) + -- ({\dx*-4.4600},{\dy*-1.0566}) + -- ({\dx*-4.4400},{\dy*-1.0496}) + -- ({\dx*-4.4200},{\dy*-1.0390}) + -- ({\dx*-4.4000},{\dy*-1.0251}) + -- ({\dx*-4.3800},{\dy*-1.0080}) + -- ({\dx*-4.3600},{\dy*-0.9878}) + -- ({\dx*-4.3400},{\dy*-0.9647}) + -- ({\dx*-4.3200},{\dy*-0.9390}) + -- ({\dx*-4.3000},{\dy*-0.9110}) + -- ({\dx*-4.2800},{\dy*-0.8810}) + -- ({\dx*-4.2600},{\dy*-0.8495}) + -- ({\dx*-4.2400},{\dy*-0.8169}) + -- ({\dx*-4.2200},{\dy*-0.7841}) + -- ({\dx*-4.2000},{\dy*-0.7518}) + -- ({\dx*-4.1800},{\dy*-0.7215}) + -- ({\dx*-4.1600},{\dy*-0.6947}) + -- ({\dx*-4.1400},{\dy*-0.6742}) + -- ({\dx*-4.1200},{\dy*-0.6644}) + -- ({\dx*-4.1000},{\dy*-0.6730}) + -- ({\dx*-4.0800},{\dy*-0.7150}) + -- ({\dx*-4.0600},{\dy*-0.8253}) + -- ({\dx*-4.0400},{\dy*-1.1085}) + -- ({\dx*-4.0200},{\dy*-2.0855}) + -- ({\dx*-4.0010},{\dy*-41.6072})} +\def\gammasinsix{({\dx*-5.9998},{\dy*6.9477}) + -- ({\dx*-5.9800},{\dy*0.1349}) + -- ({\dx*-5.9600},{\dy*0.1629}) + -- ({\dx*-5.9400},{\dy*0.2134}) + -- ({\dx*-5.9200},{\dy*0.2691}) + -- ({\dx*-5.9000},{\dy*0.3260}) + -- ({\dx*-5.8800},{\dy*0.3829}) + -- ({\dx*-5.8600},{\dy*0.4391}) + -- ({\dx*-5.8400},{\dy*0.4940}) + -- ({\dx*-5.8200},{\dy*0.5472}) + -- ({\dx*-5.8000},{\dy*0.5985}) + -- ({\dx*-5.7800},{\dy*0.6477}) + -- ({\dx*-5.7600},{\dy*0.6945}) + -- ({\dx*-5.7400},{\dy*0.7387}) + -- ({\dx*-5.7200},{\dy*0.7800}) + -- ({\dx*-5.7000},{\dy*0.8184}) + -- ({\dx*-5.6800},{\dy*0.8537}) + -- ({\dx*-5.6600},{\dy*0.8856}) + -- ({\dx*-5.6400},{\dy*0.9142}) + -- ({\dx*-5.6200},{\dy*0.9392}) + -- ({\dx*-5.6000},{\dy*0.9606}) + -- ({\dx*-5.5800},{\dy*0.9783}) + -- ({\dx*-5.5600},{\dy*0.9923}) + -- ({\dx*-5.5400},{\dy*1.0024}) + -- ({\dx*-5.5200},{\dy*1.0086}) + -- ({\dx*-5.5000},{\dy*1.0109}) + -- ({\dx*-5.4800},{\dy*1.0094}) + -- ({\dx*-5.4600},{\dy*1.0039}) + -- ({\dx*-5.4400},{\dy*0.9947}) + -- ({\dx*-5.4200},{\dy*0.9816}) + -- ({\dx*-5.4000},{\dy*0.9648}) + -- ({\dx*-5.3800},{\dy*0.9443}) + -- ({\dx*-5.3600},{\dy*0.9203}) + -- ({\dx*-5.3400},{\dy*0.8929}) + -- ({\dx*-5.3200},{\dy*0.8621}) + -- ({\dx*-5.3000},{\dy*0.8283}) + -- ({\dx*-5.2800},{\dy*0.7914}) + -- ({\dx*-5.2600},{\dy*0.7519}) + -- ({\dx*-5.2400},{\dy*0.7098}) + -- ({\dx*-5.2200},{\dy*0.6655}) + -- ({\dx*-5.2000},{\dy*0.6193}) + -- ({\dx*-5.1800},{\dy*0.5717}) + -- ({\dx*-5.1600},{\dy*0.5230}) + -- ({\dx*-5.1400},{\dy*0.4741}) + -- ({\dx*-5.1200},{\dy*0.4260}) + -- ({\dx*-5.1000},{\dy*0.3804}) + -- ({\dx*-5.0800},{\dy*0.3405}) + -- ({\dx*-5.0600},{\dy*0.3135}) + -- ({\dx*-5.0400},{\dy*0.3204}) + -- ({\dx*-5.0200},{\dy*0.4657}) + -- ({\dx*-5.0002},{\dy*41.6531})} diff --git a/buch/papers/laguerre/images/gammaplot.pdf b/buch/papers/laguerre/images/gammaplot.pdf 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/gammaplot.tex b/buch/papers/laguerre/images/gammaplot.tex new file mode 100644 index 0000000..5a68f0a --- /dev/null +++ b/buch/papers/laguerre/images/gammaplot.tex @@ -0,0 +1,73 @@ +% +% gammaplot.tex -- template for standalon tikz images +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\input{gammapaths.tex} +\begin{document} +\def\skala{1} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\definecolor{mainColor}{HTML}{D72864} % OST pink + +\draw[->] (-6.1,0) -- (5.3,0) coordinate[label={$z$}]; +\draw[->] (0,-5.1) -- (0,6.4) coordinate[label={right:$\Gamma(z)$}]; + +\foreach \x in {-1,-2,-3,-4,-5,-6}{ + \draw (\x,-0.1) -- (\x,0.1); + \draw[line width=0.1pt] (\x,-5) -- (\x,6.2); +} +\foreach \x in {1,2,3,4,5}{ + \draw (\x,-0.1) -- (\x,0.1); + \node at (\x,0) [below] {$\x$}; +} +\foreach \y in {-5,-4,-3,-2,-1,1,2,3,4,5,6}{ + \draw (-0.1,\y) -- (0.1,\y); +} +\foreach \y in {1,2,3,4,5,6}{ + \node at (0,\y) [left] {$\y$}; +} +\foreach \y in {-1,-2,-3,-4,-5}{ + \node at (0,\y) [right] {$\y$}; +} +\foreach \x in {-1,-3,-5}{ + \node at (\x,0) [below left] {$\x$}; +} +\foreach \x in {-2,-4,-6}{ + \node at (\x,0) [above left] {$\x$}; +} + +\def\dx{1} +\def\dy{1} + +\begin{scope} +\clip (-6.1,-5) rectangle (4.3,6.2); + +% \draw[color=darkgreen,line width=1.4pt] \gammasinplus; +% \draw[color=darkgreen,line width=1.4pt] \gammasinone; +% \draw[color=darkgreen,line width=1.4pt] \gammasintwo; +% \draw[color=darkgreen,line width=1.4pt] \gammasinthree; +% \draw[color=darkgreen,line width=1.4pt] \gammasinfour; +% \draw[color=darkgreen,line width=1.4pt] \gammasinfive; +% \draw[color=darkgreen,line width=1.4pt] \gammasinsix; + +\draw[color=mainColor,line width=1.4pt] \gammaplus; +\draw[color=mainColor,line width=1.4pt] \gammaone; +\draw[color=mainColor,line width=1.4pt] \gammatwo; +\draw[color=mainColor,line width=1.4pt] \gammathree; +\draw[color=mainColor,line width=1.4pt] \gammafour; +\draw[color=mainColor,line width=1.4pt] \gammafive; +\draw[color=mainColor,line width=1.4pt] \gammasix; + +\end{scope} + +\end{tikzpicture} +\end{document} + diff --git a/buch/papers/laguerre/images/integrand.pdf b/buch/papers/laguerre/images/integrand.pdf 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/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 1fe0f8b..57a6560 100644 --- a/buch/papers/laguerre/main.tex +++ b/buch/papers/laguerre/main.tex @@ -8,13 +8,32 @@ \begin{refsection} \chapterauthor{Patrik Müller} -Hier kommt eine Einleitung. +{\parindent0pt Die} Laguerre\--Polynome, +benannt nach Edmond Laguerre (1834 - 1886), +sind Lösungen der ebenfalls nach Laguerre benannten Differentialgleichung. +Laguerre entdeckte diese Polynome, als er Approximations\-methoden +für das Integral +% $\int_0^\infty \exp(-x) / x \, dx $ +\begin{align*} +\int_0^\infty \frac{e^{-x}}{x} \, dx +\end{align*} +suchte. +Darum möchten wir uns in diesem Kapitel, +ganz im Sinne des Entdeckers, +den Laguerre-Polynomen für Approximationen von Integralen mit +exponentiell-abfallenden Funktionen widmen. +Namentlich werden wir versuchen, mittels Laguerre-Polynomen und +der Gauss-Quadratur eine geeignete Approximation für die Gamma-Funktion zu finden. + +Laguerre-Polynome tauchen zudem auch in der Quantenmechanik im radialen Anteil +der Lösung für die Schrödinger-Gleichung eines Wasserstoffatoms auf. \input{papers/laguerre/definition} \input{papers/laguerre/eigenschaften} \input{papers/laguerre/quadratur} -\input{papers/laguerre/transformation} -\input{papers/laguerre/wasserstoff} +\input{papers/laguerre/gamma} +% \input{papers/laguerre/transformation} +% \input{papers/laguerre/wasserstoff} \printbibliography[heading=subbibliography] \end{refsection} diff --git a/buch/papers/laguerre/packages.tex b/buch/papers/laguerre/packages.tex index ab55228..a80d091 100644 --- a/buch/papers/laguerre/packages.tex +++ b/buch/papers/laguerre/packages.tex @@ -6,5 +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 8ab1af5..a494362 100644 --- a/buch/papers/laguerre/quadratur.tex +++ b/buch/papers/laguerre/quadratur.tex @@ -3,27 +3,185 @@ % % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % -\section{Gauss-Laguerre Quadratur -\label{laguerre:section:quadratur}} - +\section{Gauss-Quadratur + \label{laguerre:section:quadratur}} +Die Gauss-Quadratur ist ein numerisches Integrationsverfahren, +welches die Eigenschaften von orthogonalen Polynomen verwendet. +Herleitungen und Analysen der Gauss-Quadratur können im +Abschnitt~\ref{buch:orthogonal:section:gauss-quadratur} gefunden werden. +Als grundlegende Idee wird die Beobachtung, +dass viele Funktionen sich gut mit Polynomen approximieren lassen, +verwendet. +Stellt man also sicher, +dass ein Verfahren gut für Polynome funktioniert, +sollte es auch für andere Funktionen angemessene Resultate liefern. +Es wird ein Polynom verwendet, +welches an den Punkten $x_0 < x_1 < \ldots < x_n$ +die Funktionwerte~$f(x_i)$ annimmt. +Als Resultat kann das Integral via einer gewichteten Summe der Form \begin{align} - \int_a^b f(x) w(x) - \approx - \sum_{i=1}^N f(x_i) A_i - \label{laguerre:gaussquadratur} +\int_a^b f(x) w(x) \, dx +\approx +\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}} +Wir möchten nun die Gauss-Quadratur auf die Berechnung +von uneigentlichen Integralen erweitern, +spezifisch auf das Interval $(0, \infty)$. +Mit dem vorher beschriebenen Verfahren ist dies nicht direkt möglich. +Mit einer Transformation die das unendliche Intervall $(a, \infty)$ mit +\begin{align*} +x += +a + \frac{1 - t}{t} +\end{align*} +auf das Intervall $[0, 1]$ transformiert, +kann dies behoben werden. +Für unseren Fall gilt $a = 0$. +Das Integral eines Polynomes in diesem Intervall ist immer divergent. +Darum müssen wir das Polynom mit einer Funktion multiplizieren, +die schneller als jedes Polynom gegen $0$ geht, +damit das Integral immer noch konvergiert. +Die Laguerre-Polynome $L_n$ schaffen hier Abhilfe, +da ihre Gewichtsfunktion $w(x) = e^{-x}$ schneller +gegen $0$ konvergiert als jedes Polynom. +% In unserem Falle möchten wir die Gauss Quadratur auf die Laguerre-Polynome +% $L_n$ ausweiten. +% Diese sind orthogonal im Intervall $(0, \infty)$ bezüglich +% der Gewichtsfunktion $e^{-x}$. +Die Gleichung~\eqref{laguerre:gaussquadratur} lässt sich wie folgt +umformulieren: \begin{align} - \int_{0}^{\infty} f(x) e^{-x} dx - \approx - \sum_{i=1}^{N} f(x_i) A_i - \label{laguerre:laguerrequadratur} +\int_{0}^{\infty} f(x) e^{-x} dx +\approx +\sum_{i=1}^{n} f(x_i) A_i +\label{laguerre:laguerrequadratur} \end{align} +\subsubsection{Stützstellen und Gewichte} +Nach der Definition der Gauss-Quadratur müssen als Stützstellen die Nullstellen +des verwendeten Polynoms genommen werden. +Für das Laguerre-Polynom $L_n$ müssen demnach dessen Nullstellen $x_i$ und +als Gewichte $A_i$ die Integrale $l_i(x)e^{-x}$ verwendet werden. +Dabei sind +\begin{align*} +l_i(x_j) += +\delta_{ij} += +\begin{cases} +1 & i=j \\ +0 & \text{sonst} +\end{cases} +% . +\end{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 += +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} - A_i - = - \frac{x_i}{(n + 1)^2 \left[ L_{n + 1}(x_i)\right]^2} - \label{laguerre:quadratur_gewichte} +\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 += +\sum_{i=1}^n f(x_i) A_i + R_n +\end{align*} +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{laguerre:lag_error} +\end{align} +an. +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 caf270f..d21009b 100644 --- a/buch/papers/laguerre/references.bib +++ b/buch/papers/laguerre/references.bib @@ -3,33 +3,48 @@ % % (c) 2020 Autor, Hochschule Rapperswil % - -@online{laguerre:bibtex, - title = {BibTeX}, - url = {https://de.wikipedia.org/wiki/BibTeX}, - date = {2020-02-06}, - year = {2020}, - month = {2}, - day = {6} +@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{laguerre:numerical-analysis, - title = {Numerical Analysis}, - author = {David Kincaid and Ward Cheney}, - publisher = {American Mathematical Society}, - year = {2002}, - isbn = {978-8-8218-4788-6}, - inseries = {Pure and applied undegraduate texts}, - volume = {2} +@book{laguerre:abramowitz+stegun, + added-at = {2008-06-25T06:25:58.000+0200}, + address = {New York}, + author = {Abramowitz, Milton and Stegun, Irene A.}, + biburl = {https://www.bibsonomy.org/bibtex/223ec744709b3a776a1af0a3fd65cd09f/a_olympia}, + description = {BibTeX - Wikipedia, the free encyclopedia}, + edition = {ninth Dover printing, tenth GPO printing}, + interhash = {d4914a420f489f7c5129ed01ec3cf80c}, + intrahash = {23ec744709b3a776a1af0a3fd65cd09f}, + keywords = {Handbook}, + publisher = {Dover}, + pages = {890}, + timestamp = {2008-06-25T06:25:58.000+0200}, + title = {Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables}, + year = 1972 } -@article{laguerre:mendezmueller, - author = { Tabea Méndez and Andreas Müller }, - title = { Noncommutative harmonic analysis and image registration }, - journal = { Appl. Comput. Harmon. Anal.}, - year = 2019, - volume = 47, - pages = {607--627}, - url = {https://doi.org/10.1016/j.acha.2017.11.004} +@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.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_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") diff --git a/buch/papers/laguerre/transformation.tex b/buch/papers/laguerre/transformation.tex deleted file mode 100644 index 4de86b6..0000000 --- a/buch/papers/laguerre/transformation.tex +++ /dev/null @@ -1,31 +0,0 @@ -% -% transformation.tex -% -% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule -% -\section{Laguerre Transformation -\label{laguerre:section:transformation}} -\begin{align} - L \left\{ f(x) \right\} - = - \tilde{f}_\alpha(n) - = - \int_0^\infty e^{-x} x^\alpha L_n^\alpha(x) f(x) dx - \label{laguerre:transformation} -\end{align} - -\begin{align} - L^{-1} \left\{ \tilde{f}_\alpha(n) \right\} - = - f(x) - = - \sum_{n=0}^{\infty} - \begin{pmatrix} - n + \alpha \\ - n - \end{pmatrix}^{-1} - \frac{1}{\Gamma(\alpha + 1)} - \tilde{f}_\alpha(n) - L_n^\alpha(x) - \label{laguerre:inverse_transformation} -\end{align}
\ No newline at end of file diff --git a/buch/papers/laguerre/wasserstoff.tex b/buch/papers/laguerre/wasserstoff.tex deleted file mode 100644 index caaa6af..0000000 --- a/buch/papers/laguerre/wasserstoff.tex +++ /dev/null @@ -1,29 +0,0 @@ -% -% wasserstoff.tex -% -% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule -% -\section{Radialer Schwingungsanteil eines Wasserstoffatoms -\label{laguerre:section:radial_h_atom}} - -\begin{align} - \nonumber - - \frac{\hbar^2}{2m} - & - \left( - \frac{1}{r^2} \pdv{}{r} - \left( r^2 \pdv{}{r} \right) - + - \frac{1}{r^2 \sin \vartheta} \pdv{}{\vartheta} - \left( \sin \vartheta \pdv{}{\vartheta} \right) - + - \frac{1}{r^2 \sin^2 \vartheta} \pdv[2]{}{\varphi} - \right) - u(r, \vartheta, \varphi) - \\ - & - - \frac{e^2}{4 \pi \epsilon_0 r} u(r, \vartheta, \varphi) - = - E u(r, \vartheta, \varphi) - \label{laguerre:pdg_h_atom} -\end{align} |