diff options
author | Samuel Niederer <43746162+samnied@users.noreply.github.com> | 2022-07-24 12:17:00 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2022-07-24 12:17:00 +0200 |
commit | efe7c35759afb5cbae3c1683873c5159be0be09f (patch) | |
tree | 84f2e8510132352f9943bddc577ccf32cd46f2dc /buch/papers/laguerre/scripts | |
parent | add current work (diff) | |
parent | Merge pull request #26 from p1mueller/master (diff) | |
download | SeminarSpezielleFunktionen-efe7c35759afb5cbae3c1683873c5159be0be09f.tar.gz SeminarSpezielleFunktionen-efe7c35759afb5cbae3c1683873c5159be0be09f.zip |
Merge branch 'AndreasFMueller:master' into master
Diffstat (limited to 'buch/papers/laguerre/scripts')
-rw-r--r-- | buch/papers/laguerre/scripts/estimates.py | 49 | ||||
-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_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 |
11 files changed, 622 insertions, 0 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.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_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") |