diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/papers/laguerre/Makefile | 6 | ||||
| -rw-r--r-- | buch/papers/laguerre/Makefile.inc | 4 | ||||
| -rw-r--r-- | buch/papers/laguerre/definition.tex | 174 | ||||
| -rw-r--r-- | buch/papers/laguerre/eigenschaften.tex | 94 | ||||
| -rw-r--r-- | buch/papers/laguerre/images/laguerre_polynomes.pdf | bin | 0 -> 16239 bytes | |||
| -rw-r--r-- | buch/papers/laguerre/main.tex | 4 | ||||
| -rw-r--r-- | buch/papers/laguerre/packages.tex | 1 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 395 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/laguerre_plot.py | 100 | ||||
| -rw-r--r-- | buch/papers/laguerre/transformation.tex | 31 | ||||
| -rw-r--r-- | buch/papers/laguerre/wasserstoff.tex | 29 |
11 files changed, 734 insertions, 104 deletions
diff --git a/buch/papers/laguerre/Makefile b/buch/papers/laguerre/Makefile index 606d7e1..0f0985a 100644 --- a/buch/papers/laguerre/Makefile +++ b/buch/papers/laguerre/Makefile @@ -4,6 +4,8 @@ # (c) 2020 Prof Dr Andreas Mueller # -images: - @echo "no images to be created in laguerre" +images: images/laguerre_polynomes.pdf + +images/laguerre_polynomes.pdf: scripts/laguerre_plot.py + python3 scripts/laguerre_plot.py diff --git a/buch/papers/laguerre/Makefile.inc b/buch/papers/laguerre/Makefile.inc index 1eb5034..aae51f9 100644 --- a/buch/papers/laguerre/Makefile.inc +++ b/buch/papers/laguerre/Makefile.inc @@ -9,8 +9,6 @@ dependencies-laguerre = \ papers/laguerre/references.bib \ papers/laguerre/definition.tex \ papers/laguerre/eigenschaften.tex \ - papers/laguerre/quadratur.tex \ - papers/laguerre/transformation.tex \ - papers/laguerre/wasserstoff.tex + papers/laguerre/quadratur.tex diff --git a/buch/papers/laguerre/definition.tex b/buch/papers/laguerre/definition.tex index 5f6d8bd..edd2b7b 100644 --- a/buch/papers/laguerre/definition.tex +++ b/buch/papers/laguerre/definition.tex @@ -4,45 +4,151 @@ % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % \section{Definition -\label{laguerre:section:definition}} + \label{laguerre:section:definition}} \rhead{Definition} - +Die verallgemeinerte Laguerre-Differentialgleichung ist gegeben durch \begin{align} - x y''(x) + (1 - x) y'(x) + n y(x) - = - 0 - \label{laguerre:dgl} +x y''(x) + (\nu + 1 - x) y'(x) + n y(x) += +0 +, \quad +n \in \mathbb{N}_0 +, \quad +x \in \mathbb{R} +. +\label{laguerre:dgl} \end{align} - +Hier wird die verallgemeinerte Laguerre-Differentialgleichung verwendet, +weil die Lösung gleich 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 +Potenzreihenansatz. +Da wir bereits wissen, dass die Lösung orthogonale Polynome sind, +erscheint dieser Ansatz sinnvoll. +Setzt man nun den Ansatz +\begin{align*} +y(x) + & = +\sum_{k=0}^\infty a_k x^k +\\ +y'(x) + & = +\sum_{k=1}^\infty k a_k x^{k-1} += +\sum_{k=0}^\infty (k+1) a_{k+1} x^k +\\ +y''(x) + & = +\sum_{k=2}^\infty k (k-1) a_k x^{k-2} += +\sum_{k=1}^\infty (k+1) k a_{k+1} x^{k-1} +\end{align*} +in die Differentialgleichung ein, erhält man: +\begin{align*} +\sum_{k=1}^\infty (k+1) k a_{k+1} x^k ++ +(\nu + 1)\sum_{k=0}^\infty (k+1) a_{k+1} x^k +- +\sum_{k=0}^\infty k a_k x^k ++ +n \sum_{k=0}^\infty a_k x^k + & = +0 \\ +\sum_{k=1}^\infty +\left[ (k+1) k a_{k+1} + (\nu + 1)(k+1) a_{k+1} - k a_k + n a_k \right] x^k + & = +0. +\end{align*} +Daraus lässt sich die Rekursionsbeziehung +\begin{align*} +a_{k+1} + & = +\frac{k-n}{(k+1) (k + \nu + 1)} a_k +\end{align*} +ableiten. +Für ein konstantes $n$ erhalten wir als Potenzreihenlösung ein Polynom vom Grad +$n$, +denn für $k=n$ wird $a_{n+1} = 0$ und damit auch $a_{n+2}=a_{n+3}=\ldots=0$. +Aus der Rekursionsbeziehung ist zudem ersichtlich, +dass $a_0 \neq 0$ beliebig gewählt werden kann. +Wählen wir nun $a_0 = 1$, dann folgt für die Koeffizienten $a_1, a_2, a_3$ +\begin{align*} +a_1 += +-\frac{n}{1 \cdot (\nu + 1)} +, & & +a_2 += +\frac{(n-1)n}{1 \cdot 2 \cdot (\nu + 1)(\nu + 2)} +, & & +a_3 += +-\frac{(n-2)(n-1)n}{1 \cdot 2 \cdot 3 \cdot (\nu + 1)(\nu + 2)(\nu + 3)} +\end{align*} +und allgemein +\begin{align*} +k + & \leq +n: + & +a_k + & = +(-1)^k \frac{n!}{(n-k)!} \frac{1}{k!(\nu + 1)_k} += +\frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k} +\\ +k & >n: + & +a_k + & = +0. +\end{align*} +Somit erhalten wir für $\nu = 0$ die Laguerre-Polynome \begin{align} - L_n(x) - = - \sum_{k=0}^{n} - \frac{(-1)^k}{k!} - \begin{pmatrix} - n \\ - k - \end{pmatrix} - x^k - \label{laguerre:polynom} +L_n(x) += +\sum_{k=0}^{n} \frac{(-1)^k}{k!} \binom{n}{k} x^k +\label{laguerre:polynom} \end{align} - +und mit $\nu \in \mathbb{R}$ die verallgemeinerten Laguerre-Polynome \begin{align} - x y''(x) + (\alpha + 1 - x) y'(x) + n y(x) - = - 0 - \label{laguerre:generell_dgl} +L_n^\nu(x) += +\sum_{k=0}^{n} \frac{(-1)^k}{(\nu + 1)_k} \binom{n}{k} x^k. +\label{laguerre:allg_polynom} \end{align} +Durch die analytische Fortsetzung erhalten wir zudem noch die zweite Lösung der +Differentialgleichung mit der Form +\begin{align*} +\Xi_n(x) += +L_n(x) \ln(x) + \sum_{k=1}^\infty d_k x^k +\end{align*} +Nach einigen mühsamen Rechnungen, +die den Rahmen dieses Kapitel sprengen würden, +erhalten wir +\begin{align*} +\Xi_n += +L_n(x) \ln(x) ++ +\sum_{k=1}^n \frac{(-1)^k}{k!} \binom{n}{k} +(\alpha_{n-k} - \alpha_n - 2 \alpha_k)x^k ++ +(-1)^n \sum_{k=1}^\infty \frac{(k-1)!n!}{((n+k)!)^2} x^{n+k}, +\end{align*} +wobei $\alpha_0 = 0$ und $\alpha_k =\sum_{i=1}^k i^{-1}$, +$\forall k \in \mathbb{N}$. +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} -\begin{align} - L_n^\alpha (x) - = - \sum_{k=0}^{n} - \frac{(-1)^k}{k!} - \begin{pmatrix} - n + \alpha \\ - n - k - \end{pmatrix} - x^k - \label{laguerre:polynom} -\end{align} +% https://www.math.kit.edu/iana1/lehre/hm3phys2012w/media/laguerre.pdf +% http://www.physics.okayama-u.ac.jp/jeschke_homepage/E4/kapitel4.pdf diff --git a/buch/papers/laguerre/eigenschaften.tex b/buch/papers/laguerre/eigenschaften.tex index b7597e5..c589c92 100644 --- a/buch/papers/laguerre/eigenschaften.tex +++ b/buch/papers/laguerre/eigenschaften.tex @@ -4,5 +4,95 @@ % (c) 2022 Patrik Müller, Ostschweizer Fachhochschule % \section{Eigenschaften -\label{laguerre:section:eigenschaften}} -\rhead{Eigenschaften}
\ No newline at end of file + \label{laguerre:section:eigenschaften}} +\rhead{Eigenschaften} + +\subsection{Orthogonalität} +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 +Abschnitt~\ref{buch:integrale:subsection:sturm-liouville-problem}). +Der Sturm-Liouville-Operator hat die Form +\begin{align} +S += +\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right). +\label{laguerre:slop} +\end{align} +Aus der Beziehung +\begin{align} +S + & = +\Lambda +\nonumber +\\ +\frac{1}{w(x)} \left(-\frac{d}{dx}p(x) \frac{d}{dx} + q(x) \right) + & = +x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx} +\label{laguerre:sl-lag} +\end{align} +lässt sich sofort erkennen, dass $q(x) = 0$. +Ausserdem ist ersichtlich, dass $p(x)$ die Differentialgleichung +\begin{align*} +x \frac{dp}{dx} += +-(\nu + 1 - x) p, +\end{align*} +erfüllen muss. +Durch Separation erhalten wir dann +\begin{align*} +\int \frac{dp}{p} + & = +-\int \frac{\nu + 1 - x}{x}dx +\\ +\log p + & = +-\log \nu + 1 - x + C +\\ +p(x) + & = +-C x^{\nu + 1} e^{-x} +\end{align*} +Eingefügt in Gleichung~\eqref{laguerre:sl-lag} erhalten wir +\begin{align*} +\frac{C}{w(x)} +\left( +x^{\nu+1} e^{-x} \frac{d^2}{dx^2} + +(\nu + 1 - x) x^{\nu} e^{-x} \frac{d}{dx} +\right) += +x \frac{d^2}{dx^2} + (\nu + 1 - x) \frac{d}{dx}. +\end{align*} +Mittels Koeffizientenvergleich kann nun abgelesen werden, dass $w(x) = x^\nu +e^{-x}$ und $C=1$ mit $\nu > -1$. +Die Gewichtsfunktion $w(x)$ wächst für $x\rightarrow-\infty$ sehr schnell an, +deshalb ist die Laguerre-Gewichtsfunktion nur geeignet für den +Definitionsbereich $(0, \infty)$. +Bleibt nur noch sicherzustellen, dass die Randbedingungen, +\begin{align} +k_0 y(0) + h_0 p(0)y'(0) + & = +0 +\label{laguerre:sllag_randa} +\\ +k_\infty y(\infty) + h_\infty p(\infty) y'(\infty) + & = +0 +\label{laguerre:sllag_randb} +\end{align} +mit $|k_i|^2 + |h_i|^2 \neq 0,\,\forall i \in \{0, \infty\}$, erfüllt sind. +Am linken Rand (Gleichung~\eqref{laguerre:sllag_randa}) kann $y(0) = 1$, $k_0 = +0$ und $h_0 = 1$ verwendet werden, +was auch die Laguerre-Polynome ergeben haben. +Für den rechten Rand ist die Bedingung (Gleichung~\eqref{laguerre:sllag_randb}) +\begin{align*} +\lim_{x \rightarrow \infty} p(x) y'(x) + & = +\lim_{x \rightarrow \infty} -x^{\nu + 1} e^{-x} y'(x) += +0 +\end{align*} +für beliebige Polynomlösungen erfüllt für $k_\infty=0$ und $h_\infty=1$. +Damit können wir schlussfolgern, dass die Laguerre-Polynome orthogonal +bezüglich des Skalarproduktes mit der Laguerre\--Gewichtsfunktion sind. diff --git a/buch/papers/laguerre/images/laguerre_polynomes.pdf b/buch/papers/laguerre/images/laguerre_polynomes.pdf Binary files differnew file mode 100644 index 0000000..3976bc7 --- /dev/null +++ b/buch/papers/laguerre/images/laguerre_polynomes.pdf diff --git a/buch/papers/laguerre/main.tex b/buch/papers/laguerre/main.tex index 1fe0f8b..3db67d5 100644 --- a/buch/papers/laguerre/main.tex +++ b/buch/papers/laguerre/main.tex @@ -13,8 +13,8 @@ Hier kommt eine Einleitung. \input{papers/laguerre/definition} \input{papers/laguerre/eigenschaften} \input{papers/laguerre/quadratur} -\input{papers/laguerre/transformation} -\input{papers/laguerre/wasserstoff} +% \input{papers/laguerre/transformation} +% \input{papers/laguerre/wasserstoff} \printbibliography[heading=subbibliography] \end{refsection} diff --git a/buch/papers/laguerre/packages.tex b/buch/papers/laguerre/packages.tex index ab55228..4ebc172 100644 --- a/buch/papers/laguerre/packages.tex +++ b/buch/papers/laguerre/packages.tex @@ -7,4 +7,3 @@ % if your paper needs special packages, add package commands as in the % following example \usepackage{derivative} - diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb new file mode 100644 index 0000000..44f3abd --- /dev/null +++ b/buch/papers/laguerre/scripts/gamma_approx.ipynb @@ -0,0 +1,395 @@ +{ + "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/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py new file mode 100644 index 0000000..b9088d0 --- /dev/null +++ b/buch/papers/laguerre/scripts/laguerre_plot.py @@ -0,0 +1,100 @@ +#!/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/transformation.tex b/buch/papers/laguerre/transformation.tex deleted file mode 100644 index 4de86b6..0000000 --- a/buch/papers/laguerre/transformation.tex +++ /dev/null @@ -1,31 +0,0 @@ -% -% transformation.tex -% -% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule -% -\section{Laguerre Transformation -\label{laguerre:section:transformation}} -\begin{align} - L \left\{ f(x) \right\} - = - \tilde{f}_\alpha(n) - = - \int_0^\infty e^{-x} x^\alpha L_n^\alpha(x) f(x) dx - \label{laguerre:transformation} -\end{align} - -\begin{align} - L^{-1} \left\{ \tilde{f}_\alpha(n) \right\} - = - f(x) - = - \sum_{n=0}^{\infty} - \begin{pmatrix} - n + \alpha \\ - n - \end{pmatrix}^{-1} - \frac{1}{\Gamma(\alpha + 1)} - \tilde{f}_\alpha(n) - L_n^\alpha(x) - \label{laguerre:inverse_transformation} -\end{align}
\ No newline at end of file diff --git a/buch/papers/laguerre/wasserstoff.tex b/buch/papers/laguerre/wasserstoff.tex deleted file mode 100644 index caaa6af..0000000 --- a/buch/papers/laguerre/wasserstoff.tex +++ /dev/null @@ -1,29 +0,0 @@ -% -% wasserstoff.tex -% -% (c) 2022 Patrik Müller, Ostschweizer Fachhochschule -% -\section{Radialer Schwingungsanteil eines Wasserstoffatoms -\label{laguerre:section:radial_h_atom}} - -\begin{align} - \nonumber - - \frac{\hbar^2}{2m} - & - \left( - \frac{1}{r^2} \pdv{}{r} - \left( r^2 \pdv{}{r} \right) - + - \frac{1}{r^2 \sin \vartheta} \pdv{}{\vartheta} - \left( \sin \vartheta \pdv{}{\vartheta} \right) - + - \frac{1}{r^2 \sin^2 \vartheta} \pdv[2]{}{\varphi} - \right) - u(r, \vartheta, \varphi) - \\ - & - - \frac{e^2}{4 \pi \epsilon_0 r} u(r, \vartheta, \varphi) - = - E u(r, \vartheta, \varphi) - \label{laguerre:pdg_h_atom} -\end{align} |