diff options
Diffstat (limited to 'buch/papers/laguerre')
42 files changed, 3022 insertions, 616 deletions
diff --git a/buch/papers/laguerre/Makefile b/buch/papers/laguerre/Makefile index 0f0985a..85a1b83 100644 --- a/buch/papers/laguerre/Makefile +++ b/buch/papers/laguerre/Makefile @@ -3,9 +3,42 @@ # # (c) 2020 Prof Dr Andreas Mueller # +IMGFOLDER := images +PRESFOLDER := presentation -images: images/laguerre_polynomes.pdf +FIGURES := \ + images/targets.pdf \ + images/rel_error_complex.pdf \ + images/estimates.pdf \ + images/integrand.pdf \ + images/integrand_exp.pdf \ + images/laguerre_poly.pdf \ + images/rel_error_mirror.pdf \ + images/rel_error_range.pdf \ + images/rel_error_shifted.pdf \ + images/rel_error_simple.pdf \ + images/gammaplot.pdf -images/laguerre_polynomes.pdf: scripts/laguerre_plot.py - python3 scripts/laguerre_plot.py +.PHONY: all +all: images presentation +.PHONY: images +images: $(FIGURES) + +.PHONY: presentation +presentation: $(PRESFOLDER)/presentation.pdf + +images/%.pdf images/%.pgf: scripts/%.py scripts/gamma_approx.py + python3 $< + +images/gammaplot.pdf: images/gammaplot.tex images/gammapaths.tex + cd $(IMGFOLDER) && latexmk -quiet -pdf gammaplot.tex + +$(PRESFOLDER)/%.pdf: $(PRESFOLDER)/%.tex $(FIGURES) + cd $(PRESFOLDER) && latexmk -quiet -pdf $(<F) + +.PHONY: clean +clean: + rm $(FIGURES) + cd $(IMGFOLDER) && latexmk -C + cd $(PRESFOLDER) && latexmk -C
\ No newline at end of file diff --git a/buch/papers/laguerre/Makefile.inc b/buch/papers/laguerre/Makefile.inc index 12b0935..39b5d6f 100644 --- a/buch/papers/laguerre/Makefile.inc +++ b/buch/papers/laguerre/Makefile.inc @@ -10,6 +10,4 @@ dependencies-laguerre = \ papers/laguerre/definition.tex \ papers/laguerre/eigenschaften.tex \ papers/laguerre/quadratur.tex \ - papers/laguerre/gamma.tex - - + papers/laguerre/gamma.tex diff --git a/buch/papers/laguerre/definition.tex b/buch/papers/laguerre/definition.tex index d111f6f..4729a93 100644 --- a/buch/papers/laguerre/definition.tex +++ b/buch/papers/laguerre/definition.tex @@ -15,14 +15,17 @@ x y''(x) + (\nu + 1 - x) y'(x) + n y(x) 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 der selben Methode berechnet werden kann, -aber man zusätzlich die Lösung für den allgmeinen Fall erhält. -Zur Lösung der Gleichung \eqref{laguerre:dgl} verwenden wir einen +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. @@ -44,7 +47,7 @@ y''(x) = \sum_{k=1}^\infty (k+1) k a_{k+1} x^{k-1} \end{align*} -in die Differentialgleichung ein, erhält man: +in die Differentialgleichung ein, erhält man \begin{align*} \sum_{k=1}^\infty (k+1) k a_{k+1} x^k + @@ -118,6 +121,15 @@ L_n^\nu(x) \sum_{k=0}^{n} \frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k} x^k. \label{laguerre:allg_polynom} \end{align} +Die Laguerre-Polynome von Grad $0$ bis $7$ sind in +Abbildung~\ref{laguerre:fig:polyeval} dargestellt. +\begin{figure} +\centering +% \scalebox{0.8}{\input{papers/laguerre/images/laguerre_poly.pgf}} +\includegraphics[width=0.9\textwidth]{papers/laguerre/images/laguerre_poly.pdf} +\caption{Laguerre-Polynome vom Grad $0$ bis $7$} +\label{laguerre:fig:polyeval} +\end{figure} \subsection{Analytische Fortsetzung} Durch die analytische Fortsetzung erhalten wir zudem noch die zweite Lösung der @@ -126,8 +138,10 @@ Differentialgleichung mit der Form \Xi_n(x) = L_n(x) \ln(x) + \sum_{k=1}^\infty d_k x^k +. \end{align*} -Nach einigen mühsamen Rechnungen, +Nach einigen aufwändigen Rechnungen, +% die am besten ein Computeralgebrasystem übernimmt, die den Rahmen dieses Kapitel sprengen würden, erhalten wir \begin{align*} @@ -142,16 +156,5 @@ L_n(x) \ln(x) \end{align*} wobei $\alpha_0 = 0$ und $\alpha_k =\sum_{i=1}^k i^{-1}$, $\forall k \in \mathbb{N}$. -Die Laguerre-Polynome von Grad $0$ bis $7$ sind in -Abbildung~\ref{laguerre:fig:polyeval} dargestellt. -\begin{figure} -\centering -\includegraphics[width=0.7\textwidth]{% - papers/laguerre/images/laguerre_polynomes.pdf% -} -\caption{Laguerre-Polynome vom Grad $0$ bis $7$} -\label{laguerre:fig:polyeval} -\end{figure} - % https://www.math.kit.edu/iana1/lehre/hm3phys2012w/media/laguerre.pdf % http://www.physics.okayama-u.ac.jp/jeschke_homepage/E4/kapitel4.pdf diff --git a/buch/papers/laguerre/eigenschaften.tex b/buch/papers/laguerre/eigenschaften.tex index b0cc3a3..4adbe86 100644 --- a/buch/papers/laguerre/eigenschaften.tex +++ b/buch/papers/laguerre/eigenschaften.tex @@ -3,35 +3,31 @@ % % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % -\section{Eigenschaften - \label{laguerre:section:eigenschaften}} -{ -\large \color{red} -TODO: -Evtl. nur Orthogonalität hier behandeln, da nur diese für die Gauss-Quadratur -benötigt wird. -} - -Die Laguerre-Polynome besitzen einige interessante Eigenschaften -\rhead{Eigenschaften} - -\subsection{Orthogonalität - \label{laguerre:subsection:orthogonal}} -Im Abschnitt~\ref{laguerre:section:definition} haben wir behauptet, -dass die Laguerre-Polynome orthogonale Polynome sind. +\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 die Laguerre\--Differentialgleichung in ein -Sturm\--Liouville\--Problem umwandeln können, haben wir bewiesen, dass es sich -bei -den Laguerre\--Polynomen um orthogonale Polynome handelt (siehe +Wenn wir \eqref{laguerre:dgl} in ein +Sturm-Liouville-Problem umwandeln können, haben wir bewiesen, dass es sich +bei den Laguerre-Polynomen um orthogonale Polynome handelt (siehe Abschnitt~\ref{buch:integrale:subsection:sturm-liouville-problem}). -Der Sturm-Liouville-Operator hat die Form +Der Beweis kann äquivalent auch über den Sturm-Liouville-Operator \begin{align} S = \frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right). \label{laguerre:slop} \end{align} +und den Laguerre-Operator +\begin{align} +\Lambda += +x \frac{d}{dx^2} + (\nu + 1 -x) \frac{d}{dx} +\end{align} +erhalten werden, +indem wir diese Operatoren einander gleichsetzen. Aus der Beziehung \begin{align} S @@ -49,24 +45,27 @@ Ausserdem ist ersichtlich, dass $p(x)$ die Differentialgleichung \begin{align*} x \frac{dp}{dx} = --(\nu + 1 - x) p, +-(\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}{x} \, dx += +-\int \frac{\nu + 1}{x} \, dx - \int 1\, dx \\ \log p & = --\log \nu + 1 - x + C +-(\nu + 1)\log x - x + c \\ p(x) & = -C x^{\nu + 1} e^{-x} +. \end{align*} -Eingefügt in Gleichung~\eqref{laguerre:sl-lag} erhalten wir +Eingefügt in Gleichung~\eqref{laguerre:sl-lag} ergibt sich \begin{align*} \frac{C}{w(x)} \left( @@ -106,14 +105,9 @@ Für den rechten Rand ist die Bedingung (Gleichung~\eqref{laguerre:sllag_randb}) 0 \end{align*} für beliebige Polynomlösungen erfüllt für $k_\infty=0$ und $h_\infty=1$. -Damit können wir schlussfolgern, dass die Laguerre-Polynome orthogonal -bezüglich des Skalarproduktes auf dem Intervall $(0, \infty)$ mit der Laguerre\--Gewichtsfunktion -$w(x)=x^\nu e^{-x}$ sind. - - -\subsection{Rodrigues-Formel} - -\subsection{Drei-Terme Rekursion} - -\subsection{Beziehung mit der Hypergeometrischen Funktion} - +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 index e3838b0..2e5fc06 100644 --- a/buch/papers/laguerre/gamma.tex +++ b/buch/papers/laguerre/gamma.tex @@ -8,8 +8,8 @@ Die Gauss-Laguerre-Quadratur kann nun verwendet werden, um exponentiell abfallende Funktionen im Definitionsbereich $(0, \infty)$ zu berechnen. -Dabei bietet sich z.B. die Gamma-Funkion bestens an, wie wir in den folgenden -Abschnitten sehen werden. +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 @@ -19,58 +19,514 @@ Integral der Form \begin{align} \Gamma(z) & = -\int_0^\infty t^{z-1} e^{-t} dt +\int_0^\infty x^{z-1} e^{-x} \, dx , \quad \text{wobei Realteil von $z$ grösser als $0$} -, \label{laguerre:gamma} +. \end{align} -welches alle Eigenschaften erfüllt, um mit der Gauss-Laguerre-Quadratur -berechnet zu werden. +Der Term $e^{-t}$ 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 Funktionalgleichung besagt +Die Gamma-Funktion besitzt die gleiche Rekursionsbeziehung wie die Fakultät, +nämlich \begin{align} -z \Gamma(z) = \Gamma(z+1). +z \Gamma(z) += +\Gamma(z+1) +. \label{laguerre:gamma_funktional} \end{align} -Mittels dieser Gleichung kann der Wert an einer bestimmten, -geeigneten Stelle evaluiert werden und dann zurückverschoben werden, -um das gewünschte Resultat zu erhalten. + +\subsubsection{Reflektionsformel} +Die Reflektionsformel +\begin{align} +\Gamma(z) \Gamma(1 - z) += +\frac{\pi}{\sin \pi z} +,\quad +\text{für } +z \notin \mathbb{Z} +\label{laguerre:gamma_refform} +\end{align} +stellt eine Beziehung zwischen den zwei Punkten, +die aus der Spiegelung an der Geraden $\real z = 1/2$ hervorgehen, +her. +Dadurch lassen Werte der Gamma-Funktion sich für $z$ in der rechten Halbebene +leicht in die linke Halbebene übersetzen und umgekehrt. \subsection{Berechnung mittels Gauss-Laguerre-Quadratur} +In den vorherigen Abschnitten haben wir gesehen, +dass sich die Gamma-Funktion bestens für die Gauss-Laguerre-Quadratur eignet. +Nun bieten sich uns zwei Optionen, +diese zu berechnen: +\begin{enumerate} +\item Wir verwenden die verallgemeinerten Laguerre-Polynome, dann $f(x)=1$. +\item Wir verwenden die Laguerre-Polynome, dann $f(x)=x^{z-1}$. +\end{enumerate} +Die erste Variante wäre optimal auf das Problem angepasst, +allerdings müssten die Gewichte und Nullstellen für jedes $z$ +neu berechnet werden, +da sie per Definition von $z$ abhängen. +Dazu kommt, +dass die Berechnung der Gewichte $A_i$ nach \cite{laguerre:Cassity1965AbcissasCA} +\begin{align*} +A_i += +\frac{ +\Gamma(n) \Gamma(n+\nu) +} +{ +(n+\nu) +\left[L_{n-1}^{\nu}(x_i)\right]^2 +} +\end{align*} +Evaluationen der Gamma-Funktion benötigen. +Somit ist diese Methode eindeutig nicht geeignet für unser Vorhaben. + +Bei der zweiten Variante benötigen wir keine Neuberechung der Gewichte +und Nullstellen für unterschiedliche $z$. +In \eqref{laguerre:quadratur_gewichte} ist ersichtlich, +dass die Gewichte einfach zu berechnen sind. +Auch die Nullstellen können vorgängig, +mittels eines geeigneten Verfahrens, +aus den Polynomen bestimmt werden. +Als problematisch könnte sich höchstens +die zu integrierende Funktion $f(x)=x^{z-1}$ für $|z| \gg 0$ erweisen. +Somit entscheiden wir uns aufgrund der vorherigen Punkte, +die zweite Variante weiterzuverfolgen. + +\subsubsection{Direkter Ansatz} +Wenden wir also die Gauss-Laguerre-Quadratur aus +\eqref{laguerre:laguerrequadratur} auf die Gamma-Funktion +\eqref{laguerre:gamma} an, +ergibt sich +\begin{align} +\Gamma(z) +\approx +\sum_{i=1}^n x_i^{z-1} A_i. +\label{laguerre:naive_lag} +\end{align} + +\begin{figure} +\centering +% \input{papers/laguerre/images/rel_error_simple.pgf} +\includegraphics{papers/laguerre/images/rel_error_simple.pdf} +%\vspace{-12pt} +\caption{Relativer Fehler des direkten Ansatzes +für verschiedene reele Werte von $z$ und Grade $n$ der Laguerre-Polynome} +\label{laguerre:fig:rel_error_simple} +\end{figure} -Fehlerterm: +Bevor wir die Gauss-Laguerre-Quadratur anwenden, +möchten wir als ersten Schritt eine Fehlerabschätzung durchführen. +Für den Fehlerterm \eqref{laguerre:lag_error} wird die $2n$-te Ableitung +der zu integrierenden Funktion $f(\xi)$ benötigt. +Für das Integral der Gamma-Funktion ergibt sich also \begin{align*} +\frac{d^{2n}}{d\xi^{2n}} f(\xi) + & = +\frac{d^{2n}}{d\xi^{2n}} \xi^{z-1} +\\ + & = +(z - 2n)_{2n} \xi^{z - 2n - 1} +. +\end{align*} +Eingesetzt im Fehlerterm \eqref{laguerre:lag_error} resultiert +\begin{align} R_n = (z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z-2n-1} -\end{align*} +, +\label{laguerre:gamma_err_simple} +\end{align} +wobei $\xi$ ein geeigneter Wert im Interval $(0, \infty)$ ist +und $n$ der Grad des verwendeten Laguerre-Polynoms. +Eine Fehlerabschätzung mit dem Fehlerterm stellt sich als unnütz heraus, +da $R_n$ für $z < 2n - 1$ bei $\xi \rightarrow 0$ eine Singularität aufweist +und für $z > 2n - 1$ bei $\xi \rightarrow \infty$ divergiert. +Nur für den unwahrscheinlichen Fall $ z = 2n - 1$ +wäre eine Fehlerabschätzung plausibel. -\subsubsection{Finden der optimalen Berechnungsstelle} +Wenden wir nun also direkt die Gauss-Laguerre-Quadratur auf die Gamma-Funktion +an. +Dazu benötigen wir die Gewichte nach +\eqref{laguerre:quadratur_gewichte} +und als Stützstellen die Nullstellen des Laguerre-Polynomes $L_n$. +Evaluieren wir den relativen Fehler unserer Approximation zeigt sich ein +Bild wie in Abbildung~\ref{laguerre:fig:rel_error_simple}. +Man kann sehen, +wie der relative Fehler Nullstellen aufweist für ganzzahlige $z \leq 2n$. +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 an einer geeigneten Stelle evaluiert und -dann zurückverschiebt mit der Funktionalgleichung. -Dazu wollen wir den Fehlerterm in -Gleichung~\eqref{laguerre:lagurre:lag_error} anpassen und dann minimieren. -Zunächst wollen wir dies nur für $z\in \mathbb{R}$ und $0<z<1$ definieren. -Zudem nehmen wir an, dass die optimale Stelle $x^* \in \mathbb{R}$, $z < x^*$ -ist. -Dann fügen wir einen Verschiebungsterm um $m$ Stellen ein, daraus folgt +wenn man das Problem in einem geeigneten Intervall $[a(n), a(n)+1]$, +$a(n) \in \mathbb{Z}$, +evaluiert und dann mit der +Funktionalgleichung \eqref{laguerre:gamma_funktional} zurückverschiebt. +Für dieses Vorhaben führen wir einen Verschiebungsterm $m \in \mathbb{Z}$ ein. +Passen wir \eqref{laguerre:naive_lag} +mit dem Verschiebungsterm $m$ +%,der $z$ and die Stelle $z_m = z + m$ verschiebt, +an, +ergibt sich +\begin{align} +\Gamma(z) +\approx +s(z, m) \sum_{i=1}^n x_i^{z + m - 1} A_i +% && +% \text{mit } +% s(z, m) +% = +% \begin{cases} +% \displaystyle +% \frac{1}{(z - m)_m} & \text{wenn } m \geq 0\\ +% (z + m)_{-m} & \text{wenn } m < 0 +% \end{cases} +% . +\label{laguerre:shifted_lag} +\end{align} +mit \begin{align*} -R_n +s(z, m) = -\frac{(z - 2n)_{2n}}{(z - m)_m} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1} +\begin{cases} +\displaystyle +\frac{1}{(z)_m} & \text{wenn } m \geq 0 \\ +(z + m)_{-m} & \text{wenn } m < 0 +\end{cases} . \end{align*} -{ -\large \color{red} -TODO: -Geeignete Minimierung für Fehler finden, so dass sie mit den emprisich -bestimmen optimalen Punkten übereinstimmen. -} +\subsubsection{Finden der optimalen Berechnungsstelle} +Um die optimale Stelle $z^*(n) \in \left[a(n), a(n) + 1\right]$, +$z^*(n) \in \mathbb{R}$, +zu finden, +erweitern wir denn Fehlerterm \eqref{laguerre:gamma_err_simple} +und erhalten +\begin{align} +R_{n,m}(\xi) += +s(z, m) \cdot (z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z + m - 2n - 1} +,\quad +\text{für } +\xi \in (0, \infty) +\label{laguerre:gamma_err_shifted} +. +\end{align} + +\begin{figure} +\centering +\includegraphics{papers/laguerre/images/targets.pdf} +% %\vspace{-12pt} +\caption{$a$ in Abhängigkeit von $z$ und $n$} +\label{laguerre:fig:targets} +\end{figure} +% wobei ist +% mit $z^*(n) \in \mathbb{R}$ wollen wir finden, +% in dem wir den Fehlerterm \eqref{laguerre:lag_error} anpassen +% und in einem nächsten Schritt minimieren. +% Zudem nehmen wir an, +% dass $z < z^*(n)$ ist. +% Wir fügen einen Verschiebungsterm um $m \in \mathbb{N}$ Stellen ein, +% daraus folgt +% +% Damit wir den idealen Verschiebungsterm $m^*$ finden können, +% müssen wir mittels des Fehlerterms \eqref{laguerre:gamma_err_shifted} +% ein Optimierungsproblem +% +% Das Optimierungsproblem daraus lässt sich als +Daraus formulieren wir das Optimierungproblem +\begin{align*} +m^* += +\operatorname*{argmin}_m \max_\xi R_{n,m}(\xi) +. +\end{align*} +Allerdings ist die Funktion $R_{n,m}(\xi)$ unbeschränkt und +hat die gleichen Probleme wie die Fehlerabschätzung des direkten Ansatzes. +Dazu müssten wir $\xi$ versuchen, +unter Kontrolle zu bringen, +was ein äussersts schwieriges Unterfangen zu sein scheint. +Da die Gauss-Quadratur aber sowieso +nur wirklich praktisch sinnvoll für kleine $n$ ist, +können die Intervalle $[a(n), a(n)+1]$ empirisch gesucht werden. + +Wir bestimmen nun die optimalen Verschiebungsterme empirisch +für $n = 2,\ldots, 12$ im Intervall $z \in (0, 1)$, +da $z$ sowieso um den Term $m$ verschoben wird, +reicht die $m^*$ nur in diesem Intervall zu analysieren. +In Abbildung~\ref{laguerre:fig:targets} sind die empirisch bestimmten $m^*$ +abhängig von $z$ und $n$ dargestellt. +In $n$-Richtung lässt sich eine klare lineare Abhängigkeit erkennen und +die Beziehung zu $z$ ist negativ, +d.h. wenn $z$ grösser ist, wird $m^*$ kleiner. +Allerdings ist die genaue Beziehung zu $z$ +aus dieser Grafik nicht offensichtlich, +aber sie scheint regelmässig zu sein. +Es lässt die Vermutung aufkommen, +dass die Restriktion von $m^* \in \mathbb{Z}$ Rundungsprobleme verursacht. +Wir versuchen, +dieses Problem via lineare Regression und geeignete Rundung zu beheben. +Den linearen Regressor +\begin{align*} +\hat{m} += +\alpha n + \beta +\end{align*} +machen wir nur abhängig von $n$ +in dem wir den Mittelwert $\overline{m}$ von $m^*$ über $z$ berechnen. + +\begin{figure} +\centering +% \input{papers/laguerre/images/estimates.pgf} +\includegraphics{papers/laguerre/images/estimates.pdf} +%\vspace{-12pt} +\caption{Schätzung Mittelwert von $m$ und Fehler} +\label{laguerre:fig:schaetzung} +\end{figure} -\subsection{Resultate} +In Abbildung~\ref{laguerre:fig:schaetzung} sind die Resultate +der linearen Regression aufgezeigt mit $\alpha = 1.34094$ und $\beta = +0.854093$. +Die lineare Beziehung ist ganz klar ersichtlich und der Fit scheint zu genügen. +Der 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/laguerre_polynomes.pdf b/buch/papers/laguerre/images/laguerre_polynomes.pdf Binary files differdeleted file mode 100644 index 3976bc7..0000000 --- a/buch/papers/laguerre/images/laguerre_polynomes.pdf +++ /dev/null diff --git a/buch/papers/laguerre/images/rel_error_complex.pdf b/buch/papers/laguerre/images/rel_error_complex.pdf Binary files differnew file mode 100644 index 0000000..c7bb37a --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_complex.pdf diff --git a/buch/papers/laguerre/images/rel_error_mirror.pdf b/buch/papers/laguerre/images/rel_error_mirror.pdf Binary files differnew file mode 100644 index 0000000..8a27d41 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_mirror.pdf diff --git a/buch/papers/laguerre/images/rel_error_range.pdf b/buch/papers/laguerre/images/rel_error_range.pdf Binary files differnew file mode 100644 index 0000000..bb8a2d7 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_range.pdf diff --git a/buch/papers/laguerre/images/rel_error_shifted.pdf b/buch/papers/laguerre/images/rel_error_shifted.pdf Binary files differnew file mode 100644 index 0000000..b7e72dc --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_shifted.pdf diff --git a/buch/papers/laguerre/images/rel_error_simple.pdf b/buch/papers/laguerre/images/rel_error_simple.pdf Binary files differnew file mode 100644 index 0000000..3212e42 --- /dev/null +++ b/buch/papers/laguerre/images/rel_error_simple.pdf diff --git a/buch/papers/laguerre/images/targets.pdf b/buch/papers/laguerre/images/targets.pdf Binary files differnew file mode 100644 index 0000000..9514a6d --- /dev/null +++ b/buch/papers/laguerre/images/targets.pdf diff --git a/buch/papers/laguerre/main.tex b/buch/papers/laguerre/main.tex index 00e3b43..57a6560 100644 --- a/buch/papers/laguerre/main.tex +++ b/buch/papers/laguerre/main.tex @@ -8,7 +8,25 @@ \begin{refsection} \chapterauthor{Patrik Müller} -{\large \color{red} TODO: Einleitung} +{\parindent0pt Die} Laguerre\--Polynome, +benannt nach Edmond Laguerre (1834 - 1886), +sind Lösungen der ebenfalls nach Laguerre benannten Differentialgleichung. +Laguerre entdeckte diese Polynome, als er 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} diff --git a/buch/papers/laguerre/packages.tex b/buch/papers/laguerre/packages.tex index 4ebc172..a80d091 100644 --- a/buch/papers/laguerre/packages.tex +++ b/buch/papers/laguerre/packages.tex @@ -6,4 +6,4 @@ % if your paper needs special packages, add package commands as in the % following example -\usepackage{derivative} +\DeclareMathOperator{\real}{Re}
\ No newline at end of file diff --git a/buch/papers/laguerre/presentation/presentation.pdf b/buch/papers/laguerre/presentation/presentation.pdf Binary files differnew file mode 100644 index 0000000..3d00de3 --- /dev/null +++ b/buch/papers/laguerre/presentation/presentation.pdf diff --git a/buch/papers/laguerre/presentation/presentation.tex b/buch/papers/laguerre/presentation/presentation.tex new file mode 100644 index 0000000..3db69f5 --- /dev/null +++ b/buch/papers/laguerre/presentation/presentation.tex @@ -0,0 +1,134 @@ +\documentclass[ngerman, aspectratio=169, xcolor={rgb}]{beamer} + +% style +\mode<presentation>{ + \usetheme{Frankfurt} +} +%packages +\usepackage[utf8]{inputenc}\DeclareUnicodeCharacter{2212}{-} +\usepackage[ngerman]{babel} +\usepackage{graphicx} +\usepackage{array} + +\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}} +\usepackage{ragged2e} + +\usepackage{bm} % bold math +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{mathtools} +\usepackage{amsmath} +\usepackage{multirow} % multi row in tables +\usepackage{booktabs} %toprule midrule bottomrue in tables +\usepackage{scrextend} +\usepackage{textgreek} +\usepackage[rgb]{xcolor} + +\usepackage{ marvosym } % \Lightning + +\usepackage{multimedia} % embedded videos + +\usepackage{tikz} +\usepackage{pgf} +\usepackage{pgfplots} + +\usepackage{algorithmic} + +%citations +\usepackage[style=verbose,backend=biber]{biblatex} +\addbibresource{references.bib} + + +%math font +\usefonttheme[onlymath]{serif} + +%Beamer Template modifications +%\definecolor{mainColor}{HTML}{0065A3} % HSR blue +\definecolor{mainColor}{HTML}{D72864} % OST pink +\definecolor{invColor}{HTML}{28d79b} % OST pink +\definecolor{dgreen}{HTML}{38ad36} % Dark green + +%\definecolor{mainColor}{HTML}{000000} % HSR blue +\setbeamercolor{palette primary}{bg=white,fg=mainColor} +\setbeamercolor{palette secondary}{bg=orange,fg=mainColor} +\setbeamercolor{palette tertiary}{bg=yellow,fg=red} +\setbeamercolor{palette quaternary}{bg=mainColor,fg=white} %bg = Top bar, fg = active top bar topic +\setbeamercolor{structure}{fg=black} % itemize, enumerate, etc (bullet points) +\setbeamercolor{section in toc}{fg=black} % TOC sections +\setbeamertemplate{section in toc}[sections numbered] +\setbeamertemplate{subsection in toc}{% + \hspace{1.2em}{$\bullet$}~\inserttocsubsection\par} + +\setbeamertemplate{itemize items}[circle] +\setbeamertemplate{description item}[circle] +\setbeamertemplate{title page}[default][colsep=-4bp,rounded=true] +\beamertemplatenavigationsymbolsempty + +\setbeamercolor{footline}{fg=gray} +\setbeamertemplate{footline}{% + \hfill\usebeamertemplate***{navigation symbols} + \hspace{0.5cm} + \insertframenumber{}\hspace{0.2cm}\vspace{0.2cm} +} + +\usepackage{caption} +\captionsetup{labelformat=empty} + +%Title Page +\title{Laguerre-Polynome} +\subtitle{Anwendung: Approximation der Gamma-Funktion} +\author{Patrik Müller} +% \institute{OST Ostschweizer Fachhochschule} +% \institute{\includegraphics[scale=0.3]{../img/ost_logo.png}} +\date{\today} + +\input{../packages.tex} + +\newcommand*{\QED}{\hfill\ensuremath{\blacksquare}}% + +\newcommand*{\HL}{\textcolor{mainColor}} +\newcommand*{\RD}{\textcolor{red}} +\newcommand*{\BL}{\textcolor{blue}} +\newcommand*{\GN}{\textcolor{dgreen}} + +\definecolor{darkgreen}{rgb}{0,0.6,0} + + +\makeatletter +\newcount\my@repeat@count +\newcommand{\myrepeat}[2]{% + \begingroup + \my@repeat@count=\z@ + \@whilenum\my@repeat@count<#1\do{#2\advance\my@repeat@count\@ne}% + \endgroup +} +\makeatother + +\usetikzlibrary{automata,arrows,positioning,calc,shapes.geometric, fadings} + +\begin{document} + +\begin{frame} + \titlepage +\end{frame} + +\begin{frame}{Inhaltsverzeichnis} + \tableofcontents +\end{frame} + +\input{sections/laguerre} + +\input{sections/gaussquad} + +\input{sections/gamma} + +\input{sections/gamma_approx} + +\appendix +\begin{frame} + % \centering + % \Large + % \textbf{Vielen Dank für die Aufmerksamkeit} +\end{frame} + +\end{document} diff --git a/buch/papers/laguerre/presentation/sections/gamma.tex b/buch/papers/laguerre/presentation/sections/gamma.tex new file mode 100644 index 0000000..7dca39b --- /dev/null +++ b/buch/papers/laguerre/presentation/sections/gamma.tex @@ -0,0 +1,51 @@ +\section{Gamma-Funktion} + +\begin{frame}{Gamma-Funktion} +\begin{columns} + +\begin{column}{0.55\textwidth} +\begin{figure}[h] +\vspace{-16pt} +\centering +% \scalebox{0.51}{\input{../images/gammaplot.pdf}} +\includegraphics[width=1\textwidth]{../images/gammaplot.pdf} +% \caption{Gamma-Funktion} +\end{figure} +\end{column} + +\begin{column}{0.45\textwidth} +Verallgemeinerung der Fakultät +\begin{align*} +\Gamma(n) = (n-1)! +\end{align*} + +Integralformel +\begin{align*} +\Gamma(z) += +\int_0^\infty x^{z-1} e^{-x} \, dx +,\quad +\operatorname{Re} z > 0 +\end{align*} + +Funktionalgleichung +\begin{align*} +z \Gamma(z) += +\Gamma(z + 1) +\end{align*} + +Reflektionsformel +\begin{align*} +\Gamma(z) \Gamma(1 - z) += +\frac{\pi}{\sin \pi z} +, \quad +\text{für } +z \notin \mathbb{Z} +\end{align*} + +\end{column} +\end{columns} + +\end{frame}
\ No newline at end of file diff --git a/buch/papers/laguerre/presentation/sections/gamma_approx.tex b/buch/papers/laguerre/presentation/sections/gamma_approx.tex new file mode 100644 index 0000000..ecd02ab --- /dev/null +++ b/buch/papers/laguerre/presentation/sections/gamma_approx.tex @@ -0,0 +1,201 @@ +\section{Approximieren der Gamma-Funktion} + +\begin{frame}{Anwenden der Gauss-Laguerre-Quadratur auf $\Gamma(z)$} + +\begin{align*} +\Gamma(z) + & = +\int_0^\infty x^{z-1} e^{-x} \, dx +\uncover<2->{ +\approx +\sum_{i=1}^{n} f(x_i) A_i +} +\uncover<3->{ += +\sum_{i=1}^{n} x^{z-1} A_i +} +\\\\ +\uncover<4->{ + & \text{wobei } +A_i = \frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2} +\text{ und $x_i$ die Nullstellen von $L_n(x)$} +} +\end{align*} + +\end{frame} + +\begin{frame}{Fehlerabschätzung} +\begin{align*} +R_n(\xi) + & = +\frac{(n!)^2}{(2n)!} f^{(2n)}(\xi) +\\ + & = +(z - 2n)_{2n} \frac{(n!)^2}{(2n)!} \xi^{z - 2n - 1} +,\quad +0 < \xi < \infty +\end{align*} + +% \textbf{Probleme:} +\begin{itemize} +\item Funktion ist unbeschränkt +\item Maximum von $R_n$ gibt oberes Limit des Fehlers an +\uncover<2->{\item[$\Rightarrow$] Schwierig ein Maximum von $R_n(\xi)$ zu finden} +\end{itemize} +\end{frame} + +\begin{frame}{Einfacher Ansatz} + +\begin{figure}[h] +\centering +% \scalebox{0.91}{\input{../images/rel_error_simple.pgf}} +% \resizebox{!}{0.72\textheight}{\input{../images/rel_error_simple.pgf}} +\includegraphics[width=0.77\textwidth]{../images/rel_error_simple.pdf} +\caption{Relativer Fehler des einfachen Ansatzes für verschiedene reele Werte +von $z$ und Grade $n$ der Laguerre-Polynome} +\end{figure} + +\end{frame} + +\begin{frame}{Wieso sind die Resultate so schlecht?} + +\textbf{Beobachtungen} +\begin{itemize} +\item Wenn $z \in \mathbb{Z}$ relativer Fehler $\rightarrow 0$ +\item Gewisse Periodizität zu erkennen +\item Für grosse und kleine $z$ ergibt sich ein schlechter relativer Fehler +\item Es gibt Intervalle $[a,a+1]$ mit minimalem relativem Fehler +\item $a$ ist abhängig von $n$ +\end{itemize} + +\uncover<2->{ +\textbf{Ursache?} +\begin{itemize} +\item Vermutung: Integrand ist problematisch +} +\uncover<3->{ +\item[$\Rightarrow$] Analysieren von $f(x)$ und dem Integranden +} +\end{itemize} +\end{frame} + +\begin{frame}{$f(x) = x^z$} +\begin{figure}[h] +\centering +% \scalebox{0.91}{\input{../images/integrand.pgf}} +\includegraphics[width=0.8\textwidth]{../images/integrand.pdf} +% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$} +\end{figure} +\end{frame} + +\begin{frame}{Integrand $x^z e^{-x}$} +\begin{figure}[h] +\centering +% \scalebox{0.91}{\input{../images/integrand_exp.pgf}} +\includegraphics[width=0.8\textwidth]{../images/integrand_exp.pdf} +% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$} +\end{figure} +\end{frame} + +\begin{frame}{Neuer Ansatz?} + +\textbf{Vermutung} +\begin{itemize} +\item Es gibt Intervalle $[a(n), a(n)+1]$ in denen der relative Fehler minimal +ist +\item $a(n) > 0$ +\end{itemize} + +\uncover<2->{ +\textbf{Idee} +\begin{itemize} +\item[$\Rightarrow$] Berechnen von $\Gamma(z)$ im geeigneten Intervall und dann +mit Funktionalgleichung zurückverschieben +\end{itemize} +} + +\uncover<3->{ +\textbf{Wie finden wir $\boldsymbol{a(n)}$?} +\begin{itemize} +\item Minimieren des Fehlerterms mit zusätzlichem Verschiebungsterm +} +\uncover<4->{$\Rightarrow$ Schwierig das Maximum des Fehlerterms zu bestimmen} +\uncover<5->{\item Empirisch $a(n)$ bestimmen} +\uncover<6->{$\Rightarrow$ Sinnvoll, +da Gauss-Quadratur nur für kleine $n$ praktischen Nutzen hat} +\end{itemize} +\end{frame} + +\begin{frame}{Verschiebungsterm} +\begin{columns} +\begin{column}{0.625\textwidth} +\begin{figure}[h] +\centering +\includegraphics[width=1\textwidth]{../images/targets.pdf} +\caption{Optimaler Verschiebungsterm $m^*$ in Abhängigkeit von $z$ und $n$} +\end{figure} +\end{column} +\begin{column}{0.375\textwidth} +\begin{align*} +\Gamma(z) +\approx +\frac{1}{(z-m)_{m}} \sum_{i=1}^{n} x_i^{z + m - 1} A_i +\end{align*} +\end{column} +\end{columns} +\end{frame} + +\begin{frame}{Schätzen von $m^*$} +\begin{columns} +\begin{column}{0.65\textwidth} +\begin{figure} +\centering +\vspace{-12pt} +% \scalebox{0.7}{\input{../images/estimates.pgf}} +\includegraphics[width=1.0\textwidth]{../images/estimates.pdf} +% \caption{Integrand $x^z$ mit unterschiedlichen Werten für $z$} +\end{figure} +\end{column} +\begin{column}{0.34\textwidth} +\begin{align*} +\hat{m} +&= +\alpha n + \beta +\\ +&\approx +1.34093 n + 0.854093 +\\ +m^* +&= +\lceil \hat{m} - \operatorname{Re}z \rceil +\end{align*} +\end{column} +\end{columns} + +\end{frame} + +\begin{frame}{} +\begin{figure}[h] +\centering +% \scalebox{0.6}{\input{../images/rel_error_shifted.pgf}} +\includegraphics{../images/rel_error_shifted.pdf} +\caption{Relativer Fehler mit $n=8$, unterschiedlichen Verschiebungstermen $m$ und $z\in(0, 1)$} +\end{figure} +\end{frame} + +\begin{frame}{} +\begin{figure}[h] +\centering +% \scalebox{0.6}{\input{../images/rel_error_range.pgf}} +\includegraphics{../images/rel_error_range.pdf} +\caption{Relativer Fehler mit $n=8$, Verschiebungsterm $m^*$ und $z\in(-5, 5)$} +\end{figure} +\end{frame} + +\begin{frame}{Vergleich mit Lanczos-Methode} +Maximaler relativer Fehler für $n=6$ +\begin{itemize} + \item Lanczos-Methode $< 10^{-12}$ + \item Unsere Methode $\approx 10^{-6}$ +\end{itemize} +\end{frame}
\ No newline at end of file diff --git a/buch/papers/laguerre/presentation/sections/gaussquad.tex b/buch/papers/laguerre/presentation/sections/gaussquad.tex new file mode 100644 index 0000000..4d973b8 --- /dev/null +++ b/buch/papers/laguerre/presentation/sections/gaussquad.tex @@ -0,0 +1,67 @@ +\section{Gauss-Quadratur} + +\begin{frame}{Gauss-Quadratur} +\textbf{Idee} +\begin{itemize}[<+->] +\item Polynome können viele Funktionen approximieren +\item Wenn Verfahren gut für Polynome funktioniert, +sollte es auch für andere Funktionen funktionieren +\item Integrieren eines Interpolationspolynom +\item Interpolationspolynom ist durch Funktionswerte $f(x_i)$ bestimmt +$\Rightarrow$ Integral kann durch Funktionswerte berechnet werden +\item Evaluation der Funktionswerte an geeigneten Stellen +\end{itemize} +\end{frame} + +\begin{frame}{Gauss-Quadratur} +\begin{align*} +\int_{-1}^{1} f(x) \, dx +\approx +\sum_{i=1}^n f(x_i) A_i +\end{align*} + +\begin{itemize}[<+->] +\item Exakt für Polynome mit Grad $2n-1$ +\item Interpolationspolynome müssen orthogonal sein +\item Stützstellen $x_i$ sind Nullstellen des Polynoms +\item Fehler: +\begin{align*} +E += +\frac{f^{(2n)}(\xi)}{(2n)!} \int_{-1}^{1} l(x)^2 \, dx +,\quad +\text{wobei } +l(x) = \prod_{i=1}^n (x-x_i) +\end{align*} +\end{itemize} +\end{frame} + +\begin{frame}{Gauss-Laguerre-Quadratur} +\begin{itemize}[<+->] +\item Erweiterung des Integrationsintervall von $[-1, 1]$ auf $(a, b)$ +\item Hinzufügen einer Gewichtsfunktion +\item Bei uneigentlichen Integralen muss Gewichtsfunktion schneller als jedes +Integrationspolynom gegen $0$ gehen +\item[$\Rightarrow$] Für Laguerre-Polynome haben wir den Definitionsbereich +$(0, \infty)$ und die Gewichtsfunktion $w(x) = e^{-x}$ +\begin{align*} +\int_0^\infty & f(x) e^{-x} \, dx +\approx +\sum_{i=1}^n f(x_i) A_i +\\ + & \text{wobei } +A_i = \frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2} +\text{ und $x_i$ die Nullstellen von $L_n(x)$} +\end{align*} +\end{itemize} +\end{frame} + +\begin{frame}{Fehler der Gauss-Laguerre-Quadratur} +\begin{align*} +R_n += +\frac{(n!)^2}{(2n)!} f^{(2n)}(\xi) +,\quad +0 < \xi < \infty +\end{align*} +\end{frame}
\ No newline at end of file diff --git a/buch/papers/laguerre/presentation/sections/laguerre.tex b/buch/papers/laguerre/presentation/sections/laguerre.tex new file mode 100644 index 0000000..f99214e --- /dev/null +++ b/buch/papers/laguerre/presentation/sections/laguerre.tex @@ -0,0 +1,91 @@ +\section{Laguerre-Polynome} + +\begin{frame}{Laguerre-Differentialgleichung} + +\begin{itemize} +\item Benannt nach Edmond Nicolas Laguerre (1834-1886) +\item Aus Artikel von 1879, +in dem er $\int_0^\infty \exp(-x)/x \, dx$ analysierte +\end{itemize} + +\begin{align*} +x y''(x) + (1 - x) y'(x) + n y(x) + & = +0 +, \quad +n \in \mathbb{N}_0 +, \quad +x \in \mathbb{R} +\end{align*} + +\end{frame} + +\begin{frame}{Lösen der Differentialgleichung} + +\begin{align*} +x y''(x) + (1 - x) y'(x) + n y(x) + & = +0 +\\ +\end{align*} + +\uncover<2->{ +\centering +\begin{tikzpicture}[remember picture,overlay] +%% use here too +\path[draw=mainColor, very thick,->](0, 1.1) to +node[anchor=west]{Potenzreihenansatz} (0, -0.8); +\end{tikzpicture} +} + +\begin{align*} +\uncover<3->{ +L_n(x) + & = +\sum_{k=0}^{n} \frac{(-1)^k}{k!} \binom{n}{k} x^k +} +\end{align*} +\uncover<4->{ +\begin{itemize} + \item Die Lösungen der DGL sind die Laguerre-Polynome +\end{itemize} +} +\end{frame} + +\begin{frame} +\begin{figure}[h] +\centering +% \resizebox{0.74\textwidth}{!}{\input{../images/laguerre_poly.pgf}} +\includegraphics[width=0.7\textwidth]{../images/laguerre_poly.pdf} +\caption{Laguerre-Polynome vom Grad $0$ bis $7$} +\end{figure} +\end{frame} + +\begin{frame}{Orthogonalität} +\begin{itemize}[<+->] +\item Beweis: Umformen in Sturm-Liouville-Problem (siehe Paper) +\begin{alignat*}{5} +((p(x) &y'(x)))' + q(x) &y(x) +&= +\lambda &w(x) &y(x) +\\ +((x e^{-x} &y'(x)))' + 0 &y(x) +&= +n &e^{-x} &y(x) +\end{alignat*} +\item Definitionsbereich $(0, \infty)$ +\item Gewichtsfunktion $w(x) = e^{-x}$ +\end{itemize} + +\uncover<4->{ +\begin{align*} +\int_0^\infty e^{-x} L_n(x) L_m(x) \, dx += +0 +,\quad +n \neq m +,\quad +n, m \in \mathbb{N} +\end{align*} +} +\end{frame}
\ No newline at end of file diff --git a/buch/papers/laguerre/quadratur.tex b/buch/papers/laguerre/quadratur.tex index 60fad7f..a494362 100644 --- a/buch/papers/laguerre/quadratur.tex +++ b/buch/papers/laguerre/quadratur.tex @@ -5,35 +5,70 @@ % \section{Gauss-Quadratur \label{laguerre:section:quadratur}} - {\large \color{red} TODO: Einleitung und kurze Beschreibung Gauss-Quadratur} +Die Gauss-Quadratur ist ein numerisches Integrationsverfahren, +welches die Eigenschaften von orthogonalen Polynomen 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) +\int_a^b f(x) w(x) \, dx \approx -\sum_{i=1}^N f(x_i) A_i +\sum_{i=1}^n f(x_i) A_i \label{laguerre:gaussquadratur} \end{align} +berechnet werden. +Die Gauss-Quadratur ist exakt für Polynome mit Grad $2n -1$, +wenn ein Interpolationspolynom von Grad $n$ gewählt wurde. \subsection{Gauss-Laguerre-Quadratur \label{laguerre:subsection:gausslag-quadratur}} -Die Gauss-Quadratur kann auch auf Skalarprodukte mit Gewichtsfunktionen -ausgeweitet werden. -In unserem Falle möchten wir die Gauss Quadratur auf die Laguerre-Polynome -$L_n$ ausweiten. -Diese sind orthogonal im Intervall $(0, \infty)$ bezüglich -der Gewichtsfunktion $e^{-x}$. -Gleichung~\eqref{laguerre:laguerrequadratur} lässt sich wiefolgt umformulieren: +Wir möchten nun die Gauss-Quadratur auf die Berechnung +von uneigentlichen Integralen erweitern, +spezifisch auf das Interval $(0, \infty)$. +Mit dem vorher beschriebenen Verfahren ist dies nicht direkt möglich. +Mit einer Transformation die das unendliche Intervall $(a, \infty)$ mit +\begin{align*} +x += +a + \frac{1 - t}{t} +\end{align*} +auf das Intervall $[0, 1]$ transformiert, +kann dies behoben werden. +Für unseren Fall gilt $a = 0$. +Das Integral eines Polynomes in diesem Intervall ist immer divergent. +Darum müssen wir das 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 +\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. -Das heisst für das Laguerre-Polynom $L_n$ müssen dessen Nullstellen $x_i$ und -als Gewichte $A_i$ werden die Integrale $l_i(x)e^{-x}$ verwendet werden. +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) @@ -41,39 +76,112 @@ l_i(x_j) \delta_{ij} = \begin{cases} -1 & i=j \\ -0 & \text{sonst.} +1 & i=j \\ +0 & \text{sonst} \end{cases} +% . \end{align*} -Laut \cite{abramowitz+stegun} sind die Gewichte also -\begin{align} +die Lagrangschen Interpolationspolynome. +Laut \cite{laguerre:hildebrand2013introduction} können die Gewichte mit +\begin{align*} A_i + & = +-\frac{C_{n+1} \gamma_n}{C_n \phi'_n(x_i) \phi_{n+1} (x_i)} +\end{align*} +berechnet werden. +$C_i$ entspricht dabei dem Koeffizienten von $x^i$ +des orthogonalen Polynoms $\phi_n(x)$, $\forall i =0,\ldots,n$ und +\begin{align*} +\gamma_n += +\int_0^\infty w(x) \phi_n^2(x)\,dx +\end{align*} +dem Normalisierungsfaktor. +Wir setzen nun $\phi_n(x) = L_n(x)$ und +nutzen den Vorzeichenwechsel der Laguerre-Koeffizienten aus, +damit erhalten wir +\begin{align*} +A_i + & = +-\frac{C_{n+1} \gamma_n}{C_n L'_n(x_i) L_{n+1} (x_i)} +\\ + & = \frac{C_n}{C_{n-1}} \frac{\gamma_{n-1}}{L_{n-1}(x_i) L'_n(x_i)} +. +\end{align*} +Für Laguerre-Polynome gilt +\begin{align*} +\frac{C_n}{C_{n-1}} += +-\frac{1}{n} +\quad \text{und} \quad +\gamma_n = -\frac{x_i}{(n + 1)^2 \left[ L_{n + 1}(x_i)\right]^2} +1 +. +\end{align*} +Daraus folgt +\begin{align} +A_i +&= +- \frac{1}{n L_{n-1}(x_i) L'_n(x_i)} +. +\label{laguerre:gewichte_lag_temp} +\end{align} +Nun kann die Rekursionseigenschaft der Laguerre-Polynome +\begin{align*} +x L'_n(x) +&= +n L_n(x) - n L_{n-1}(x) +\\ +&= (x - n - 1) L_n(x) + (n + 1) L_{n+1}(x) +\end{align*} +umgeformt werden und da $x_i$ die Nullstellen von $L_n(x)$ sind, +vereinfacht sich der Term zu +\begin{align*} +x_i L'_n(x_i) +&= +- n L_{n-1}(x_i) +\\ +&= + (n + 1) L_{n+1}(x_i) +. +\end{align*} +Setzen wir das nun in \eqref{laguerre:gewichte_lag_temp} ein, +ergibt sich +\begin{align} +\nonumber +A_i +&= +\frac{1}{x_i \left[ L'_n(x_i) \right]^2} +\\ +&= +\frac{x_i}{(n+1)^2 \left[ L_{n+1}(x_i) \right]^2} . \label{laguerre:quadratur_gewichte} \end{align} \subsubsection{Fehlerterm} +Die Gauss-Laguerre-Quadratur mit $n$ Stützstellen berechnet Integrale +von Polynomen bis zum Grad $2n - 1$ exakt. +Für beliebige Funktionen kann eine Fehlerabschätzung angegeben werden. Der Fehlerterm $R_n$ folgt direkt aus der Approximation \begin{align*} -\int_0^{\infty} f(x) e^{-x} dx +\int_0^{\infty} f(x) e^{-x} \, dx = \sum_{i=1}^n f(x_i) A_i + R_n \end{align*} -un \cite{abramowitz+stegun} gibt in als +und \cite{laguerre:abramowitz+stegun} gibt ihn als \begin{align} R_n -= + & = +\frac{f^{(2n)}(\xi)}{(2n)!} \int_0^\infty l(x)^2 e^{-x}\,dx +\\ + & = \frac{(n!)^2}{(2n)!} f^{(2n)}(\xi) ,\quad 0 < \xi < \infty -\label{lagurre:lag_error} +\label{laguerre:lag_error} \end{align} an. - -{ -\large \color{red} -TODO: -Noch mehr Text / bessere Beschreibungen in allen Abschnitten -} +Der Fehler ist also abhängig von der $2n$-ten Ableitung +der zu integrierenden Funktion. diff --git a/buch/papers/laguerre/references.bib b/buch/papers/laguerre/references.bib index 6956ade..d21009b 100644 --- a/buch/papers/laguerre/references.bib +++ b/buch/papers/laguerre/references.bib @@ -3,8 +3,17 @@ % % (c) 2020 Autor, Hochschule Rapperswil % +@book{laguerre:hildebrand2013introduction, + title={Introduction to Numerical Analysis: Second Edition}, + author={Hildebrand, F.B.}, + isbn={9780486318554}, + series={Dover Books on Mathematics}, + year={2013}, + publisher={Dover Publications}, + pages = {389} +} -@book{abramowitz+stegun, +@book{laguerre:abramowitz+stegun, added-at = {2008-06-25T06:25:58.000+0200}, address = {New York}, author = {Abramowitz, Milton and Stegun, Irene A.}, @@ -19,4 +28,23 @@ timestamp = {2008-06-25T06:25:58.000+0200}, title = {Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables}, year = 1972 +} + +@article{laguerre:Cassity1965AbcissasCA, + title={Abcissas, coefficients, and error term for the generalized Gauss-Laguerre quadrature formula using the zero ordinate}, + author={C. Ronald Cassity}, + journal={Mathematics of Computation}, + year={1965}, + volume={19}, + pages={287-296} +} + +@online{laguerre:lanczos, + title = {Lanczos Approximation}, + author={Eric W. Weisstein}, + url = {https://mathworld.wolfram.com/LanczosApproximation.html}, + date = {2022-07-18}, + year = {2022}, + month = {7}, + day = {18} }
\ No newline at end of file diff --git a/buch/papers/laguerre/scripts/estimates.py b/buch/papers/laguerre/scripts/estimates.py new file mode 100644 index 0000000..21551f3 --- /dev/null +++ b/buch/papers/laguerre/scripts/estimates.py @@ -0,0 +1,49 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + import gamma_approx as ga + import targets + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + N = 200 + ns = np.arange(2, 13) + step = 1 / (N - 1) + x = np.linspace(step, 1 - step, N + 1) + + bests = targets.find_best_loc(N, ns=ns) + mean_m = np.mean(bests, -1) + + intercept, bias = np.polyfit(ns, mean_m, 1) + fig, axs = plt.subplots( + 2, num=1, sharex=True, clear=True, constrained_layout=True, figsize=(4.5, 3.6) + ) + xl = np.array([ns[0] - 0.5, ns[-1] + 0.5]) + axs[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$") + axs[0].plot(ns, mean_m, "x", label=r"$\overline{m}$") + axs[1].plot( + ns, ((intercept * ns + bias) - mean_m), "-x", label=r"$\hat{m} - \overline{m}$" + ) + axs[0].set_xlim(*xl) + axs[0].set_xticks(ns) + axs[0].set_yticks(np.arange(np.floor(mean_m[0]), np.ceil(mean_m[-1]) + 0.1, 2)) + # axs[0].set_title("Schätzung von Mittelwert") + # axs[1].set_title("Fehler") + axs[-1].set_xlabel(r"$n$") + for ax in axs: + ax.grid(1) + ax.legend() + # fig.savefig(f"{ga.img_path}/estimates.pgf") + fig.savefig(f"{ga.img_path}/estimates.pdf") + + print(f"Intercept={intercept:.6g}, Bias={bias:.6g}") + predicts = np.ceil(intercept * ns[:, None] + bias - np.real(x)) + print(f"Error: {np.mean(np.abs(bests - predicts))}") diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb deleted file mode 100644 index 44f3abd..0000000 --- a/buch/papers/laguerre/scripts/gamma_approx.ipynb +++ /dev/null @@ -1,395 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Gauss-Laguerre Quadratur für die Gamma-Funktion\n", - "\n", - "$$\n", - " \\Gamma(z)\n", - " = \n", - " \\int_0^\\infty t^{z-1}e^{-t}dt\n", - "$$\n", - "\n", - "$$\n", - " \\int_0^\\infty f(x) e^{-x} dx \n", - " \\approx \n", - " \\sum_{i=1}^{N} f(x_i) w_i\n", - " \\qquad\\text{ wobei }\n", - " w_i = \\frac{x_i}{(n+1)^2 [L_{n+1}(x_i)]^2}\n", - "$$\n", - "und $x_i$ sind Nullstellen des Laguerre Polynoms $L_n(x)$\n", - "\n", - "Der Fehler ist gegeben als\n", - "\n", - "$$\n", - " E \n", - " =\n", - " \\frac{(n!)^2}{(2n)!} f^{(2n)}(\\xi) \n", - " = \n", - " \\frac{(-2n + z)_{2n}}{(z-m)_m} \\frac{(n!)^2}{(2n)!} \\xi^{z + m - 2n - 1}\n", - "$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from cmath import exp, pi, sin, sqrt\n", - "import scipy.special\n", - "\n", - "EPSILON = 1e-07\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lanczos_p = [\n", - " 676.5203681218851,\n", - " -1259.1392167224028,\n", - " 771.32342877765313,\n", - " -176.61502916214059,\n", - " 12.507343278686905,\n", - " -0.13857109526572012,\n", - " 9.9843695780195716e-6,\n", - " 1.5056327351493116e-7,\n", - "]\n", - "\n", - "\n", - "def drop_imag(z):\n", - " if abs(z.imag) <= EPSILON:\n", - " z = z.real\n", - " return z\n", - "\n", - "\n", - "def lanczos_gamma(z):\n", - " z = complex(z)\n", - " if z.real < 0.5:\n", - " y = pi / (sin(pi * z) * lanczos_gamma(1 - z)) # Reflection formula\n", - " else:\n", - " z -= 1\n", - " x = 0.99999999999980993\n", - " for (i, pval) in enumerate(lanczos_p):\n", - " x += pval / (z + i + 1)\n", - " t = z + len(lanczos_p) - 0.5\n", - " y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x\n", - " return drop_imag(y)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "zeros, weights = np.polynomial.laguerre.laggauss(8)\n", - "# zeros = np.array(\n", - "# [\n", - "# 1.70279632305101000e-1,\n", - "# 9.03701776799379912e-1,\n", - "# 2.25108662986613069e0,\n", - "# 4.26670017028765879e0,\n", - "# 7.04590540239346570e0,\n", - "# 1.07585160101809952e1,\n", - "# 1.57406786412780046e1,\n", - "# 2.28631317368892641e1,\n", - "# ]\n", - "# )\n", - "\n", - "# weights = np.array(\n", - "# [\n", - "# 3.69188589341637530e-1,\n", - "# 4.18786780814342956e-1,\n", - "# 1.75794986637171806e-1,\n", - "# 3.33434922612156515e-2,\n", - "# 2.79453623522567252e-3,\n", - "# 9.07650877335821310e-5,\n", - "# 8.48574671627253154e-7,\n", - "# 1.04800117487151038e-9,\n", - "# ]\n", - "# )\n", - "\n", - "\n", - "def pochhammer(z, n):\n", - " return np.prod(z + np.arange(n))\n", - "\n", - "\n", - "def find_shift(z, target):\n", - " factor = 1.0\n", - " steps = int(np.floor(target - np.real(z)))\n", - " zs = z + steps\n", - " if steps > 0:\n", - " factor = 1 / pochhammer(z, steps)\n", - " elif steps < 0:\n", - " factor = pochhammer(zs, -steps)\n", - " return zs, factor\n", - "\n", - "\n", - "def laguerre_gamma(z, x, w, target=11):\n", - " # res = 0.0\n", - " z = complex(z)\n", - " if z.real < 1e-3:\n", - " res = pi / (\n", - " sin(pi * z) * laguerre_gamma(1 - z, x, w, target)\n", - " ) # Reflection formula\n", - " else:\n", - " z_shifted, correction_factor = find_shift(z, target)\n", - " res = np.sum(x ** (z_shifted - 1) * w)\n", - " res *= correction_factor\n", - " res = drop_imag(res)\n", - " return res\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def eval_laguerre(x, target=12):\n", - " return np.array([laguerre_gamma(xi, zeros, weights, target) for xi in x])\n", - "\n", - "\n", - "def eval_lanczos(x):\n", - " return np.array([lanczos_gamma(xi) for xi in x])\n", - "\n", - "\n", - "def eval_mean_laguerre(x, targets):\n", - " return np.mean([eval_laguerre(x, target) for target in targets], 0)\n", - "\n", - "\n", - "def calc_rel_error(x, y):\n", - " return (y - x) / x\n", - "\n", - "\n", - "def evaluate(x, target=12):\n", - " lanczos_gammas = eval_lanczos(x)\n", - " laguerre_gammas = eval_laguerre(x, target)\n", - " rel_error = calc_rel_error(lanczos_gammas, laguerre_gammas)\n", - " return rel_error\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Test with real values" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Empirische Tests zeigen:\n", - "- $n=4 \\Rightarrow m=6$\n", - "- $n=5 \\Rightarrow m=7$ oder $m=8$\n", - "- $n=6 \\Rightarrow m=9$\n", - "- $n=7 \\Rightarrow m=10$\n", - "- $n=8 \\Rightarrow m=11$ oder $m=12$\n", - "- $n=9 \\Rightarrow m=13$\n", - "- $n=10 \\Rightarrow m=14$\n", - "- $n=11 \\Rightarrow m=15$ oder $m=16$\n", - "- $n=12 \\Rightarrow m=17$\n", - "- $n=13 \\Rightarrow m=18 \\Rightarrow $ Beginnt numerisch instabil zu werden \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "zeros, weights = np.polynomial.laguerre.laggauss(12)\n", - "targets = np.arange(16, 21)\n", - "mean_targets = ((16, 17),)\n", - "x = np.linspace(EPSILON, 1 - EPSILON, 101)\n", - "_, axs = plt.subplots(\n", - " 2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n", - ")\n", - "\n", - "lanczos = eval_lanczos(x)\n", - "for mean_target in mean_targets:\n", - " vals = eval_mean_laguerre(x, mean_target)\n", - " rel_error_mean = calc_rel_error(lanczos, vals)\n", - " axs[0].plot(x, rel_error_mean, label=mean_target)\n", - " axs[1].semilogy(x, np.abs(rel_error_mean), label=mean_target)\n", - "\n", - "mins = []\n", - "maxs = []\n", - "for target in targets:\n", - " rel_error = evaluate(x, target)\n", - " mins.append(np.min(np.abs(rel_error[(0.1 <= x) & (x <= 0.9)])))\n", - " maxs.append(np.max(np.abs(rel_error)))\n", - " axs[0].plot(x, rel_error, label=target)\n", - " axs[1].semilogy(x, np.abs(rel_error), label=target)\n", - "# axs[0].set_ylim(*(np.array([-1, 1]) * 3.5e-8))\n", - "\n", - "axs[0].set_xlim(x[0], x[-1])\n", - "axs[1].set_ylim(np.min(mins), 1.04*np.max(maxs))\n", - "for ax in axs:\n", - " ax.legend()\n", - " ax.grid(which=\"both\")\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "targets = (16, 17)\n", - "xmax = 15\n", - "x = np.linspace(-xmax + EPSILON, xmax - EPSILON, 1000)\n", - "\n", - "mean_lag = eval_mean_laguerre(x, targets)\n", - "lanczos = eval_lanczos(x)\n", - "rel_error = calc_rel_error(lanczos, mean_lag)\n", - "rel_error_simple = evaluate(x, targets[-1])\n", - "# rel_error = evaluate(x, target)\n", - "\n", - "_, axs = plt.subplots(\n", - " 2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n", - ")\n", - "axs[0].plot(x, rel_error, label=targets)\n", - "axs[1].semilogy(x, np.abs(rel_error), label=targets)\n", - "axs[0].plot(x, rel_error_simple, label=targets[-1])\n", - "axs[1].semilogy(x, np.abs(rel_error_simple), label=targets[-1])\n", - "axs[0].set_xlim(x[0], x[-1])\n", - "# axs[0].set_ylim(*(np.array([-1, 1]) * 4.2e-8))\n", - "# axs[1].set_ylim(1e-10, 5e-8)\n", - "for ax in axs:\n", - " ax.legend()\n", - "\n", - "x2 = np.linspace(-5 + EPSILON, 5, 4001)\n", - "_, ax = plt.subplots(constrained_layout=True, figsize=(8, 6))\n", - "ax.plot(x2, eval_mean_laguerre(x2, targets))\n", - "ax.set_xlim(x2[0], x2[-1])\n", - "ax.set_ylim(-7.5, 25)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Test with complex values" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "targets = (16, 17)\n", - "vals = np.linspace(-5 + EPSILON, 5, 100)\n", - "x, y = np.meshgrid(vals, vals)\n", - "mesh = x + 1j * y\n", - "input = mesh.flatten()\n", - "\n", - "mean_lag = eval_mean_laguerre(input, targets).reshape(mesh.shape)\n", - "lanczos = eval_lanczos(input).reshape(mesh.shape)\n", - "rel_error = np.abs(calc_rel_error(lanczos, mean_lag))\n", - "\n", - "lag = eval_laguerre(input, targets[-1]).reshape(mesh.shape)\n", - "rel_error_simple = np.abs(calc_rel_error(lanczos, lag))\n", - "# rel_error = evaluate(x, target)\n", - "\n", - "fig, axs = plt.subplots(\n", - " 2,\n", - " 2,\n", - " sharex=True,\n", - " sharey=True,\n", - " clear=True,\n", - " constrained_layout=True,\n", - " figsize=(12, 10),\n", - ")\n", - "_c = axs[0, 1].pcolormesh(x, y, np.log10(np.abs(lanczos - mean_lag)), shading=\"gouraud\")\n", - "_c = axs[0, 0].pcolormesh(x, y, np.log10(np.abs(lanczos - lag)), shading=\"gouraud\")\n", - "fig.colorbar(_c, ax=axs[0, :])\n", - "_c = axs[1, 1].pcolormesh(x, y, np.log10(rel_error), shading=\"gouraud\")\n", - "_c = axs[1, 0].pcolormesh(x, y, np.log10(rel_error_simple), shading=\"gouraud\")\n", - "fig.colorbar(_c, ax=axs[1, :])\n", - "_ = axs[0, 0].set_title(\"Absolute Error\")\n", - "_ = axs[1, 0].set_title(\"Relative Error\")\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "z = 0.5\n", - "ns = [4, 5, 5, 6, 7, 8, 8, 9, 10, 11, 11, 12] # np.arange(4, 13)\n", - "ms = np.arange(6, 18)\n", - "xi = np.logspace(0, 2, 201)[:, None]\n", - "lanczos = eval_lanczos([z])[0]\n", - "\n", - "_, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(12, 8))\n", - "ax.grid(1)\n", - "for n, m in zip(ns, ms):\n", - " zeros, weights = np.polynomial.laguerre.laggauss(n)\n", - " c = scipy.special.factorial(n) ** 2 / scipy.special.factorial(2 * n)\n", - " e = np.abs(\n", - " scipy.special.poch(z - 2 * n, 2 * n)\n", - " / scipy.special.poch(z - m, m)\n", - " * c\n", - " * xi ** (z - 2 * n + m - 1)\n", - " )\n", - " ez = np.sum(\n", - " scipy.special.poch(z - 2 * n, 2 * n)\n", - " / scipy.special.poch(z - m, m)\n", - " * c\n", - " * zeros[:, None] ** (z - 2 * n + m - 1),\n", - " 0,\n", - " )\n", - " lag = eval_laguerre([z], m)[0]\n", - " err = np.abs(lanczos - lag)\n", - " # print(m+z,ez)\n", - " # for zi,ezi in zip(z[0], ez):\n", - " # print(f\"{m+zi}: {ezi}\")\n", - " # ax.semilogy(xi, e, color=color)\n", - " lines = ax.loglog(xi, e, label=str(n))\n", - " ax.axhline(err, color=lines[0].get_color())\n", - " # ax.set_xticks(np.arange(xi[-1] + 1))\n", - " # ax.set_ylim(1e-8, 1e5)\n", - "_ = ax.legend()\n", - "# _ = ax.legend([f\"z={zi}\" for zi in z[0]])\n", - "# _ = [ax.axvline(x) for x in zeros]\n" - ] - } - ], - "metadata": { - "interpreter": { - "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90" - }, - "kernelspec": { - "display_name": "Python 3.8.10 64-bit", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.10" - }, - "orig_nbformat": 4 - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py new file mode 100644 index 0000000..5b09e59 --- /dev/null +++ b/buch/papers/laguerre/scripts/gamma_approx.py @@ -0,0 +1,116 @@ +from pathlib import Path + +import numpy as np +import scipy.special + +EPSILON = 1e-7 +root = str(Path(__file__).parent) +img_path = f"{root}/../images" +fontsize = "medium" +cmap = "plasma" + + +def _prep_zeros_and_weights(x, w, n): + if x is None or w is None: + return np.polynomial.laguerre.laggauss(n) + return x, w + + +def drop_imag(z): + if abs(z.imag) <= EPSILON: + z = z.real + return z + + +def pochhammer(z, n): + return np.prod(z + np.arange(n)) + + +def find_optimal_shift(z, n): + mhat = 1.34093 * n + 0.854093 + steps = int(np.floor(mhat - np.real(z))) + return steps + + +def get_shifting_factor(z, steps): + factor = 1.0 + if steps > 0: + factor = 1 / pochhammer(z, steps) + elif steps < 0: + factor = pochhammer(z + steps, -steps) + return factor + + +def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11): + x, w = _prep_zeros_and_weights(x, w, n) + n = len(x) + + z += 0j + steps = int(np.floor(target - np.real(z))) + z_shifted = z + steps + correction_factor = get_shifting_factor(z, steps) + + res = np.sum(x ** (z_shifted - 1) * w) + res *= correction_factor + res = drop_imag(res) + return res + + +def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty + x, w = _prep_zeros_and_weights(x, w, n) + n = len(x) + + z += 0j + opt_shift = find_optimal_shift(z, n) + correction_factor = get_shifting_factor(z, opt_shift) + z_shifted = z + opt_shift + + res = np.sum(x ** (z_shifted - 1) * w) + res *= correction_factor + res = drop_imag(res) + return res + + +def laguerre_gamma_simple(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty + x, w = _prep_zeros_and_weights(x, w, n) + z += 0j + res = np.sum(x ** (z - 1) * w) + res = drop_imag(res) + return res + + +def laguerre_gamma_mirror(z, x=None, w=None, n=8): + if z == 0.0: + return np.infty + x, w = _prep_zeros_and_weights(x, w, n) + z += 0j + if z.real < 1e-3: + return np.pi / ( + np.sin(np.pi * z) * laguerre_gamma_simple(1 - z, x, w) + ) # Reflection formula + return laguerre_gamma_simple(z, x, w) + + +def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs): + x, w = _prep_zeros_and_weights(x, w, n) + if func == "simple": + f = laguerre_gamma_simple + elif func == "mirror": + f = laguerre_gamma_mirror + elif func == "optimal_shifted": + f = laguerre_gamma_opt_shifted + else: + f = laguerre_gamma_shifted + return np.array([f(zi, x, w, n, **kwargs) for zi in z]) + + +def calc_rel_error(x, y): + mask = np.abs(x) != np.infty + rel_error = np.zeros_like(y) + rel_error[mask] = (y[mask] - x[mask]) / x[mask] + rel_error[~mask] = 0.0 + return rel_error diff --git a/buch/papers/laguerre/scripts/integrand.py b/buch/papers/laguerre/scripts/integrand.py new file mode 100644 index 0000000..e970721 --- /dev/null +++ b/buch/papers/laguerre/scripts/integrand.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +"""Plot for integrand of gamma function with shifting terms.""" + +if __name__ == "__main__": + import os + from pathlib import Path + + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + EPSILON = 1e-12 + xlims = np.array([-3, 3]) + + root = str(Path(__file__).parent) + img_path = f"{root}/../images" + os.makedirs(img_path, exist_ok=True) + + t = np.logspace(*xlims, 1001)[:, None] + + z = np.array([-4.5, -2, -1, -0.5, 0.0, 0.5, 1, 2, 4.5]) + r = t ** z + + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(4, 2.4)) + ax.semilogx(t, r) + ax.set_xlim(*(10.0 ** xlims)) + ax.set_ylim(1e-3, 40) + ax.set_xlabel(r"$x$") + # ax.set_ylabel(r"$x^z$") + ax.grid(1, "both") + labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)] + ax.legend(labels, ncol=2, loc="upper left", fontsize="small") + fig.savefig(f"{img_path}/integrand.pdf") diff --git a/buch/papers/laguerre/scripts/integrand_exp.py b/buch/papers/laguerre/scripts/integrand_exp.py new file mode 100644 index 0000000..e649b26 --- /dev/null +++ b/buch/papers/laguerre/scripts/integrand_exp.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +"""Plot for integrand of gamma function with shifting terms.""" + +if __name__ == "__main__": + import os + from pathlib import Path + + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + EPSILON = 1e-12 + xlims = np.array([-3, 3]) + + root = str(Path(__file__).parent) + img_path = f"{root}/../images" + os.makedirs(img_path, exist_ok=True) + + t = np.logspace(*xlims, 1001)[:, None] + + z = np.array([-1, -0.5, 0.0, 0.5, 1, 2, 3, 4, 4.5]) + e = np.exp(-t) + r = t ** z * e + + fig, ax = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(4, 2.4)) + ax.semilogx(t, r) + # ax.plot(t,np.exp(-t)) + ax.set_xlim(10 ** (-2), 20) + ax.set_ylim(1e-3, 10) + ax.set_xlabel(r"$x$") + # ax.set_ylabel(r"$x^z e^{-x}$") + ax.grid(1, "both") + labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)] + ax.legend(labels, ncol=2, loc="upper left", fontsize="small") + # fig.savefig(f"{img_path}/integrand_exp.pgf") + fig.savefig(f"{img_path}/integrand_exp.pdf") + # plt.show() diff --git a/buch/papers/laguerre/scripts/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py deleted file mode 100644 index b9088d0..0000000 --- a/buch/papers/laguerre/scripts/laguerre_plot.py +++ /dev/null @@ -1,100 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding:utf-8 -*- -"""Some plots for Laguerre Polynomials.""" - -import os -from pathlib import Path - -import matplotlib.pyplot as plt -import numpy as np -import scipy.special as ss - - -def get_ticks(start, end, step=1): - ticks = np.arange(start, end, step) - return ticks[ticks != 0] - - -N = 1000 -step = 5 -t = np.linspace(-1.05, 10.5, N)[:, None] -root = str(Path(__file__).parent) -img_path = f"{root}/../images" -os.makedirs(img_path, exist_ok=True) - - -# fig = plt.figure(num=1, clear=True, tight_layout=True, figsize=(5.5, 3.7)) -# ax = fig.add_subplot(axes_class=AxesZero) -fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(6, 4)) -for n in np.arange(0, 8): - k = np.arange(0, n + 1)[None] - L = np.sum((-1) ** k * ss.binom(n, k) / ss.factorial(k) * t ** k, -1) - ax.plot(t, L, label=f"n={n}") - -ax.set_xticks(get_ticks(int(t[0]), t[-1]), minor=True) -ax.set_xticks(get_ticks(0, t[-1], step)) -ax.set_xlim(t[0], t[-1] + 0.1 * (t[1] - t[0])) -ax.set_xlabel(r"$x$", x=1.0, labelpad=-10, rotation=0, fontsize="large") - -ylim = 13 -ax.set_yticks(np.arange(-ylim, ylim), minor=True) -ax.set_yticks(np.arange(-step * (ylim // step), ylim, step)) -ax.set_ylim(-ylim, ylim) -ax.set_ylabel(r"$y$", y=0.95, labelpad=-18, rotation=0, fontsize="large") - -ax.legend(ncol=2, loc=(0.125, 0.01), fontsize="large") - -# set the x-spine -ax.spines[["left", "bottom"]].set_position("zero") -ax.spines[["right", "top"]].set_visible(False) -ax.xaxis.set_ticks_position("bottom") -hlx = 0.4 -dx = t[-1, 0] - t[0, 0] -dy = 2 * ylim -hly = dy / dx * hlx -dps = fig.dpi_scale_trans.inverted() -bbox = ax.get_window_extent().transformed(dps) -width, height = bbox.width, bbox.height - -# manual arrowhead width and length -hw = 1.0 / 60.0 * dy -hl = 1.0 / 30.0 * dx -lw = 0.5 # axis line width -ohg = 0.0 # arrow overhang - -# compute matching arrowhead length and width -yhw = hw / dy * dx * height / width -yhl = hl / dx * dy * width / height - -# draw x and y axis -ax.arrow( - t[-1, 0] - hl, - 0, - hl, - 0.0, - fc="k", - ec="k", - lw=lw, - head_width=hw, - head_length=hl, - overhang=ohg, - length_includes_head=True, - clip_on=False, -) - -ax.arrow( - 0, - ylim - yhl, - 0.0, - yhl, - fc="k", - ec="k", - lw=lw, - head_width=yhw, - head_length=yhl, - overhang=ohg, - length_includes_head=True, - clip_on=False, -) - -fig.savefig(f"{img_path}/laguerre_polynomes.pdf") diff --git a/buch/papers/laguerre/scripts/laguerre_poly.py b/buch/papers/laguerre/scripts/laguerre_poly.py new file mode 100644 index 0000000..9700ab4 --- /dev/null +++ b/buch/papers/laguerre/scripts/laguerre_poly.py @@ -0,0 +1,108 @@ +import numpy as np + + +def get_ticks(start, end, step=1): + ticks = np.arange(start, end, step) + return ticks[ticks != 0] + + +if __name__ == "__main__": + import os + from pathlib import Path + + import matplotlib as mpl + import matplotlib.pyplot as plt + import scipy.special as ss + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + N = 1000 + step = 5 + t = np.linspace(-1.05, 10.5, N)[:, None] + root = str(Path(__file__).parent) + img_path = f"{root}/../images" + os.makedirs(img_path, exist_ok=True) + + # fig = plt.figure(num=1, clear=True, tight_layout=True, figsize=(5.5, 3.7)) + # ax = fig.add_subplot(axes_class=AxesZero) + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(6, 4)) + for n in np.arange(0, 8): + k = np.arange(0, n + 1)[None] + L = np.sum((-1) ** k * ss.binom(n, k) / ss.factorial(k) * t ** k, -1) + ax.plot(t, L, label=f"$n={n}$") + + ax.set_xticks(get_ticks(int(t[0]), t[-1]), minor=True) + ax.set_xticks(get_ticks(0, t[-1], step)) + ax.set_xlim(t[0], t[-1] + 0.1 * (t[1] - t[0])) + ax.set_xlabel(r"$x$", x=1.0, labelpad=-10, rotation=0, fontsize="large") + + ylim = 13 + ax.set_yticks(get_ticks(-ylim, ylim), minor=True) + ax.set_yticks(get_ticks(-step * (ylim // step), ylim, step)) + ax.set_ylim(-ylim, ylim) + ax.set_ylabel(r"$y$", y=0.95, labelpad=-18, rotation=0, fontsize="large") + + ax.legend(ncol=2, loc=(0.125, 0.01), fontsize="large") + + # set the x-spine + ax.spines[["left", "bottom"]].set_position("zero") + ax.spines[["right", "top"]].set_visible(False) + ax.xaxis.set_ticks_position("bottom") + hlx = 0.4 + dx = t[-1, 0] - t[0, 0] + dy = 2 * ylim + hly = dy / dx * hlx + dps = fig.dpi_scale_trans.inverted() + bbox = ax.get_window_extent().transformed(dps) + width, height = bbox.width, bbox.height + + # manual arrowhead width and length + hw = 1.0 / 60.0 * dy + hl = 1.0 / 30.0 * dx + lw = 0.5 # axis line width + ohg = 0.0 # arrow overhang + + # compute matching arrowhead length and width + yhw = hw / dy * dx * height / width + yhl = hl / dx * dy * width / height + + # draw x and y axis + ax.arrow( + t[-1, 0] - hl, + 0, + hl, + 0.0, + fc="k", + ec="k", + lw=lw, + head_width=hw, + head_length=hl, + overhang=ohg, + length_includes_head=True, + clip_on=False, + ) + + ax.arrow( + 0, + ylim - yhl, + 0.0, + yhl, + fc="k", + ec="k", + lw=lw, + head_width=yhw, + head_length=yhl, + overhang=ohg, + length_includes_head=True, + clip_on=False, + ) + + # fig.savefig(f"{img_path}/laguerre_poly.pgf") + fig.savefig(f"{img_path}/laguerre_poly.pdf") + # plt.show() diff --git a/buch/papers/laguerre/scripts/rel_error_complex.py b/buch/papers/laguerre/scripts/rel_error_complex.py new file mode 100644 index 0000000..4a714fa --- /dev/null +++ b/buch/papers/laguerre/scripts/rel_error_complex.py @@ -0,0 +1,43 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + import scipy.special + + import gamma_approx as ga + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + xmax = 4 + vals = np.linspace(-xmax + ga.EPSILON, xmax, 100) + x, y = np.meshgrid(vals, vals) + mesh = x + 1j * y + input = mesh.flatten() + + lanczos = scipy.special.gamma(mesh) + lag = ga.eval_laguerre_gamma(input, n=8, func="optimal_shifted").reshape(mesh.shape) + rel_error = np.abs(ga.calc_rel_error(lanczos, lag)) + + fig, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(3.5, 2.1)) + _c = ax.pcolormesh( + x, y, rel_error, shading="gouraud", cmap=ga.cmap, norm=mpl.colors.LogNorm() + ) + cbr = fig.colorbar(_c, ax=ax) + cbr.minorticks_off() + # ax.set_title("Relative Error") + ax.set_xlabel("Re($z$)") + ax.set_ylabel("Im($z$)") + minor_ticks = np.arange(-xmax, xmax + ga.EPSILON) + ticks = np.arange(-xmax, xmax + ga.EPSILON, 2) + ax.set_xticks(ticks) + ax.set_xticks(minor_ticks, minor=True) + ax.set_yticks(ticks) + ax.set_yticks(minor_ticks, minor=True) + fig.savefig(f"{ga.img_path}/rel_error_complex.pdf") + # plt.show() diff --git a/buch/papers/laguerre/scripts/rel_error_mirror.py b/buch/papers/laguerre/scripts/rel_error_mirror.py new file mode 100644 index 0000000..7348d5e --- /dev/null +++ b/buch/papers/laguerre/scripts/rel_error_mirror.py @@ -0,0 +1,38 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + import scipy.special + + import gamma_approx as ga + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + xmin = -15 + xmax = 15 + ns = np.arange(2, 12, 2) + ylim = np.array([-11, 1]) + x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400) + gamma = scipy.special.gamma(x) + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5)) + for n in ns: + gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="mirror") + rel_err = ga.calc_rel_error(gamma, gamma_lag) + ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$") + ax.set_xlim(x[0], x[-1]) + ax.set_ylim(*(10.0 ** ylim)) + ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5)) + ax.set_xticks(np.arange(xmin, xmax), minor=True) + ax.set_yticks(10.0 ** np.arange(*ylim, 2)) + ax.set_xlabel(r"$z$") + # ax.set_ylabel("Relativer Fehler") + ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize) + ax.grid(1, "both") + # fig.savefig(f"{ga.img_path}/rel_error_mirror.pgf") + fig.savefig(f"{ga.img_path}/rel_error_mirror.pdf") diff --git a/buch/papers/laguerre/scripts/rel_error_range.py b/buch/papers/laguerre/scripts/rel_error_range.py new file mode 100644 index 0000000..ece3b6d --- /dev/null +++ b/buch/papers/laguerre/scripts/rel_error_range.py @@ -0,0 +1,41 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + import scipy.special + + import gamma_approx as ga + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + N = 1201 + xmax = 6 + xmin = -xmax + ns = np.arange(2, 12, 2) + ylim = np.array([-11, -1.2]) + + x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, N) + gamma = scipy.special.gamma(x) + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2)) + for n in ns: + gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted") + rel_err = ga.calc_rel_error(gamma, gamma_lag) + ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$") + ax.set_xlim(x[0], x[-1]) + ax.set_ylim(*(10.0 ** ylim)) + ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 2)) + ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON), minor=True) + ax.set_yticks(10.0 ** np.arange(*ylim, 2)) + ax.set_yticks(10.0 ** np.arange(*ylim, 1), "", minor=True) + ax.set_xlabel(r"$z$") + # ax.set_ylabel("Relativer Fehler") + ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize) + ax.grid(1, "both") + # fig.savefig(f"{ga.img_path}/rel_error_range.pgf") + fig.savefig(f"{ga.img_path}/rel_error_range.pdf") + # plt.show() diff --git a/buch/papers/laguerre/scripts/rel_error_shifted.py b/buch/papers/laguerre/scripts/rel_error_shifted.py new file mode 100644 index 0000000..f53c89b --- /dev/null +++ b/buch/papers/laguerre/scripts/rel_error_shifted.py @@ -0,0 +1,40 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + import scipy.special + + import gamma_approx as ga + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + n = 8 # order of Laguerre polynomial + N = 200 # number of points in interval + + step = 1 / (N - 1) + x = np.linspace(step, 1 - step, N + 1) + targets = np.arange(10, 14) + gamma = scipy.special.gamma(x) + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.1)) + for target in targets: + gamma_lag = ga.eval_laguerre_gamma(x, target=target, n=n, func="shifted") + rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lag)) + ax.semilogy(x, rel_error, label=f"$m={target}$", linewidth=3) + gamma_lgo = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted") + rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lgo)) + ax.semilogy(x, rel_error, "m", linestyle=":", label="$m^*$", linewidth=3) + ax.set_xlim(x[0], x[-1]) + ax.set_ylim(5e-9, 5e-8) + ax.set_xlabel(r"$z$") + ax.set_xticks(np.linspace(0, 1, 6)) + ax.set_xticks(np.linspace(0, 1, 11), minor=True) + ax.grid(1, "both") + ax.legend(ncol=1, fontsize=ga.fontsize) + # fig.savefig(f"{ga.img_path}/rel_error_shifted.pgf") + fig.savefig(f"{ga.img_path}/rel_error_shifted.pdf") + # plt.show() diff --git a/buch/papers/laguerre/scripts/rel_error_simple.py b/buch/papers/laguerre/scripts/rel_error_simple.py new file mode 100644 index 0000000..686500b --- /dev/null +++ b/buch/papers/laguerre/scripts/rel_error_simple.py @@ -0,0 +1,41 @@ +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + import scipy.special + + import gamma_approx as ga + + # mpl.rc("text", usetex=True) + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + # mpl.rcParams.update({"font.family": "serif", "font.serif": "TeX Gyre Termes"}) + + # Simple / naive + xmin = -5 + xmax = 30 + ns = np.arange(2, 12, 2) + ylim = np.array([-11, 6]) + x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400) + gamma = scipy.special.gamma(x) + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5)) + for n in ns: + gamma_lag = ga.eval_laguerre_gamma(x, n=n) + rel_err = ga.calc_rel_error(gamma, gamma_lag) + ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$") + ax.set_xlim(x[0], x[-1]) + ax.set_ylim(*(10.0 ** ylim)) + ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5)) + ax.set_xticks(np.arange(xmin, xmax), minor=True) + ax.set_yticks(10.0 ** np.arange(*ylim, 2)) + ax.set_xlabel(r"$z$") + # ax.set_ylabel("Relativer Fehler") + ax.legend(ncol=3, fontsize=ga.fontsize) + ax.grid(1, "both") + # fig.savefig(f"{ga.img_path}/rel_error_simple.pgf") + fig.savefig(f"{ga.img_path}/rel_error_simple.pdf") diff --git a/buch/papers/laguerre/scripts/targets.py b/buch/papers/laguerre/scripts/targets.py new file mode 100644 index 0000000..3bc7f52 --- /dev/null +++ b/buch/papers/laguerre/scripts/targets.py @@ -0,0 +1,58 @@ +import numpy as np +import scipy.special + +import gamma_approx as ga + + +def find_best_loc(N=200, a=1.375, b=0.5, ns=None): + if ns is None: + ns = np.arange(2, 13) + bests = [] + step = 1 / (N - 1) + x = np.linspace(step, 1 - step, N + 1) + gamma = scipy.special.gamma(x) + for n in ns: + zeros, weights = np.polynomial.laguerre.laggauss(n) + est = np.ceil(b + a * n) + targets = np.arange(max(est - 2, 0), est + 3) + rel_error = [] + for target in targets: + gamma_lag = ga.eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted") + rel_error.append(np.abs(ga.calc_rel_error(gamma, gamma_lag))) + rel_error = np.stack(rel_error, -1) + best = np.argmin(rel_error, -1) + targets[0] + bests.append(best) + return np.stack(bests, 0) + + +if __name__ == "__main__": + import matplotlib as mpl + import matplotlib.pyplot as plt + + mpl.rcParams.update( + { + "mathtext.fontset": "stix", + "font.family": "serif", + "font.serif": "TeX Gyre Termes", + } + ) + + N = 200 + ns = np.arange(2, 13) + + bests = find_best_loc(N, ns=ns) + + fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(3.5, 2.1)) + v = ax.imshow(bests, cmap=ga.cmap, aspect="auto", interpolation="nearest") + plt.colorbar(v, ax=ax, label=r"$m^*$") + ticks = np.arange(0, N + 1, N // 5) + ax.set_xlim(0, 1) + ax.set_xticks(ticks) + ax.set_xticklabels([f"{v:.2f}" for v in ticks / N]) + ax.set_xticks(np.arange(0, N + 1, N // 20), minor=True) + ax.set_yticks(np.arange(len(ns))) + ax.set_yticklabels(ns) + ax.set_xlabel(r"$z$") + ax.set_ylabel(r"$n$") + # fig.savefig(f"{ga.img_path}/targets.pgf") + fig.savefig(f"{ga.img_path}/targets.pdf") |