diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 80 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.py | 197 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/integrand.py | 49 | ||||
| -rw-r--r-- | buch/papers/laguerre/scripts/laguerre_plot.py | 5 |
4 files changed, 317 insertions, 14 deletions
diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb index 44f3abd..a8280aa 100644 --- a/buch/papers/laguerre/scripts/gamma_approx.ipynb +++ b/buch/papers/laguerre/scripts/gamma_approx.ipynb @@ -136,14 +136,17 @@ "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", + " # 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", + " 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" ] @@ -207,9 +210,9 @@ "metadata": {}, "outputs": [], "source": [ - "zeros, weights = np.polynomial.laguerre.laggauss(12)\n", - "targets = np.arange(16, 21)\n", - "mean_targets = ((16, 17),)\n", + "zeros, weights = np.polynomial.laguerre.laggauss(8)\n", + "targets = np.arange(9, 14)\n", + "mean_targets = ((9, 10),)\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", @@ -226,7 +229,7 @@ "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", + " mins.append(np.min(np.abs(rel_error[(0.05 <= x) & (x <= 0.95)])))\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", @@ -365,6 +368,59 @@ "# _ = ax.legend([f\"z={zi}\" for zi in z[0]])\n", "# _ = [ax.axvline(x) for x in zeros]\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bests = []\n", + "N = 200\n", + "step = 1 / (N - 1)\n", + "a = 11 / 8\n", + "b = 1 / 2\n", + "x = np.linspace(step, 1 - step, N + 1)\n", + "ns = np.arange(2, 13)\n", + "for n in ns:\n", + " zeros, weights = np.polynomial.laguerre.laggauss(n)\n", + " est = np.ceil(b + a * n)\n", + " targets = np.arange(max(est - 2, 0), est + 3)\n", + " rel_errors = np.stack([np.abs(evaluate(x, target)) for target in targets], -1)\n", + " best = np.argmin(rel_errors, -1) + targets[0]\n", + " bests.append(best)\n", + "bests = np.stack(bests, 0)\n", + "\n", + "fig, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(5, 3))\n", + "v = ax.imshow(bests, cmap=\"inferno\", aspect=\"auto\")\n", + "plt.colorbar(v, ax=ax, label=r'$m$')\n", + "ticks = np.arange(0, N + 1, 10)\n", + "ax.set_xlim(0, 1)\n", + "ax.set_xticks(ticks, [f\"{v:.2f}\" for v in ticks / N])\n", + "ax.set_xticks(np.arange(N + 1), minor=True)\n", + "ax.set_yticks(np.arange(len(ns)), ns)\n", + "ax.set_xlabel(r\"$z$\")\n", + "ax.set_ylabel(r\"$n$\")\n", + "# for best in bests:\n", + "# print(\", \".join([f\"{int(b):2d}\" for b in best]))\n", + "# print(np.unique(bests, return_counts=True))\n", + "\n", + "targets = np.mean(bests, -1)\n", + "intercept, bias = np.polyfit(ns, targets, 1)\n", + "_, axs2 = plt.subplots(2, sharex=True, clear=True, constrained_layout=True)\n", + "xl = np.array([1, ns[-1] + 1])\n", + "axs2[0].plot(ns, intercept * ns + bias)\n", + "axs2[0].plot(ns, targets, \"x\")\n", + "axs2[1].plot(ns, ((intercept * ns + bias) - targets), \"-x\")\n", + "print(np.mean(bests, -1))\n", + "print(f\"Intercept={intercept:.6g}, Bias={bias:.6g}\")\n", + "\n", + "\n", + "predicts = np.ceil(intercept * ns[:, None] + bias - x)\n", + "print(np.sum(np.abs(bests-predicts)))\n", + "# for best in predicts:\n", + "# print(\", \".join([f\"{int(b):2d}\" for b in best]))\n" + ] } ], "metadata": { diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py new file mode 100644 index 0000000..90843b1 --- /dev/null +++ b/buch/papers/laguerre/scripts/gamma_approx.py @@ -0,0 +1,197 @@ +from pathlib import Path + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import scipy.special + +EPSILON = 1e-7 +root = str(Path(__file__).parent) +img_path = f"{root}/../images" + + +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_shift(z, target): + factor = 1.0 + steps = int(np.floor(target - np.real(z))) + zs = z + steps + if steps > 0: + factor = 1 / pochhammer(z, steps) + elif steps < 0: + factor = pochhammer(zs, -steps) + return zs, factor + + +def laguerre_gamma_shift(z, x=None, w=None, n=8, target=11): + x, w = _prep_zeros_and_weights(x, w, n) + + z += 0j + z_shifted, correction_factor = find_shift(z, target) + 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): + 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): + 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 + else: + f = laguerre_gamma_shift + return np.array([f(zi, x, w, n, **kwargs) for zi in z]) + + +def calc_rel_error(x, y): + return (y - x) / x + + +ns = np.arange(2, 12, 2) + +# Simple / naive +xmin = -5 +xmax = 30 +ylim = np.array([-11, 6]) +x = np.linspace(xmin + EPSILON, xmax - 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 = eval_laguerre_gamma(x, n=n) + rel_err = 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 + EPSILON, 5)) +ax.set_xticks(np.arange(xmin, xmax), minor=True) +ax.set_yticks(10.0 ** np.arange(*ylim, 2)) +ax.set_yticks(10.0 ** np.arange(*ylim, 2)) +ax.set_xlabel(r"$z$") +ax.set_ylabel("Relativer Fehler") +ax.legend(ncol=3, fontsize="small") +ax.grid(1, "both") +fig.savefig(f"{img_path}/rel_error_simple.pgf") + + +# Mirrored +xmin = -15 +xmax = 15 +ylim = np.array([-11, 1]) +x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400) +gamma = scipy.special.gamma(x) +fig2, ax2 = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(5, 2.5)) +for n in ns: + gamma_lag = eval_laguerre_gamma(x, n=n, func="mirror") + rel_err = calc_rel_error(gamma, gamma_lag) + ax2.semilogy(x, np.abs(rel_err), label=f"$n={n}$") +ax2.set_xlim(x[0], x[-1]) +ax2.set_ylim(*(10.0 ** ylim)) +ax2.set_xticks(np.arange(xmin, xmax + EPSILON, 5)) +ax2.set_xticks(np.arange(xmin, xmax), minor=True) +ax2.set_yticks(10.0 ** np.arange(*ylim, 2)) +# locmin = mpl.ticker.LogLocator(base=10.0,subs=0.1*np.arange(1,10),numticks=100) +# ax2.yaxis.set_minor_locator(locmin) +# ax2.yaxis.set_minor_formatter(mpl.ticker.NullFormatter()) +ax2.set_xlabel(r"$z$") +ax2.set_ylabel("Relativer Fehler") +ax2.legend(ncol=1, loc="upper left", fontsize="small") +ax2.grid(1, "both") +fig2.savefig(f"{img_path}/rel_error_mirror.pgf") + + +# Move to target +bests = [] +N = 200 +step = 1 / (N - 1) +a = 11 / 8 +b = 1 / 2 +x = np.linspace(step, 1 - step, N + 1) +gamma = scipy.special.gamma(x)[:, None] +ns = np.arange(2, 13) +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) + gamma_lag = np.stack( + [ + eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted") + for target in targets + ], + -1, + ) + rel_error = np.abs(calc_rel_error(gamma, gamma_lag)) + best = np.argmin(rel_error, -1) + targets[0] + bests.append(best) +bests = np.stack(bests, 0) + +fig3, ax3 = plt.subplots(num=3, clear=True, constrained_layout=True, figsize=(5, 3)) +v = ax3.imshow(bests, cmap="inferno", aspect="auto", interpolation="nearest") +plt.colorbar(v, ax=ax3, label=r"$m$") +ticks = np.arange(0, N + 1, N // 5) +ax3.set_xlim(0, 1) +ax3.set_xticks(ticks, [f"{v:.2f}" for v in ticks / N]) +ax3.set_xticks(np.arange(0, N + 1, N // 20), minor=True) +ax3.set_yticks(np.arange(len(ns)), ns) +ax3.set_xlabel(r"$z$") +ax3.set_ylabel(r"$n$") +fig3.savefig(f"{img_path}/targets.pdf") + +targets = np.mean(bests, -1) +intercept, bias = np.polyfit(ns, targets, 1) +fig4, axs4 = plt.subplots( + 2, num=4, sharex=True, clear=True, constrained_layout=True, figsize=(5, 4) +) +xl = np.array([ns[0] - 0.5, ns[-1] + 0.5]) +axs4[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$") +axs4[0].plot(ns, targets, "x", label=r"$\bar{m}$") +axs4[1].plot( + ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \bar{m}$" +) +axs4[0].set_xlim(*xl) +# axs4[0].set_title("Schätzung von Mittelwert") +# axs4[1].set_title("Fehler") +axs4[-1].set_xlabel(r"$z$") +for ax in axs4: + ax.grid(1) + ax.legend() +fig4.savefig(f"{img_path}/schaetzung.pgf") + +print(f"Intercept={intercept:.6g}, Bias={bias:.6g}") +predicts = np.ceil(intercept * ns[:, None] + bias - x) +print(f"Error: {int(np.sum(np.abs(bests-predicts)))}") + +# plt.show() diff --git a/buch/papers/laguerre/scripts/integrand.py b/buch/papers/laguerre/scripts/integrand.py new file mode 100644 index 0000000..0cf43d1 --- /dev/null +++ b/buch/papers/laguerre/scripts/integrand.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +"""Plot for integrand of gamma function with shifting terms.""" + +import os +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +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=(5, 3)) +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}/integrands.pgf") + +z2 = np.array([-1, -0.5, 0.0, 0.5, 1, 2, 3, 4, 4.5]) +e = np.exp(-t) +r2 = t ** z2 * e + +fig2, ax2 = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(5, 3)) +ax2.semilogx(t, r2) +# ax2.plot(t,np.exp(-t)) +ax2.set_xlim(10 ** (-2), 20) +ax2.set_ylim(1e-3, 10) +ax2.set_xlabel(r"$x$") +ax2.set_ylabel(r"$x^z e^{-x}$") +ax2.grid(1, "both") +labels =[f"$z={zi: 3.1f}$" for zi in np.squeeze(z2)] +ax2.legend(labels, ncol=2, loc="upper left", fontsize="small") +fig2.savefig(f"{img_path}/integrands_exp.pgf") +# plt.show() diff --git a/buch/papers/laguerre/scripts/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py index b9088d0..1be3552 100644 --- a/buch/papers/laguerre/scripts/laguerre_plot.py +++ b/buch/papers/laguerre/scripts/laguerre_plot.py @@ -29,7 +29,7 @@ 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.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)) @@ -97,4 +97,5 @@ ax.arrow( clip_on=False, ) -fig.savefig(f"{img_path}/laguerre_polynomes.pdf") +fig.savefig(f"{img_path}/laguerre_polynomes.pgf") +# plt.show() |