diff options
author | Andreas Müller <andreas.mueller@ost.ch> | 2022-07-19 18:48:09 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2022-07-19 18:48:09 +0200 |
commit | 97448d520134eed27c72cd60221910d5d2191ec9 (patch) | |
tree | 414029c38e3faa222e9ba08645d4305996b4fd48 /buch/papers/laguerre/scripts/gamma_approx.ipynb | |
parent | makefile fix (diff) | |
parent | Correct typos, improve grammar (diff) | |
download | SeminarSpezielleFunktionen-97448d520134eed27c72cd60221910d5d2191ec9.tar.gz SeminarSpezielleFunktionen-97448d520134eed27c72cd60221910d5d2191ec9.zip |
Merge pull request #23 from p1mueller/master
Erster Entwurf Laguerre
Diffstat (limited to 'buch/papers/laguerre/scripts/gamma_approx.ipynb')
-rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 395 |
1 files changed, 0 insertions, 395 deletions
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 -} |