diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 395 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/laguerre_plot.py | 100 |
2 files changed, 495 insertions, 0 deletions
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") |