diff options
Diffstat (limited to '')
-rw-r--r-- | buch/papers/laguerre/scripts/estimates.py | 49 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 395 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.py | 116 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/integrand.py | 42 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/integrand_exp.py | 46 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/laguerre_plot.py | 100 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/laguerre_poly.py | 108 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/rel_error_complex.py | 43 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/rel_error_mirror.py | 38 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/rel_error_range.py | 41 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/rel_error_shifted.py | 40 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/rel_error_simple.py | 41 | ||||
-rw-r--r-- | buch/papers/laguerre/scripts/targets.py | 58 |
13 files changed, 622 insertions, 495 deletions
diff --git a/buch/papers/laguerre/scripts/estimates.py b/buch/papers/laguerre/scripts/estimates.py new file mode 100644 index 0000000..1acd7f7 --- /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(1, 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..05db5d3 --- /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=-14, 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..e1ea36a --- /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 = 25 + 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..69f94ba --- /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(1, 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") |