aboutsummaryrefslogtreecommitdiffstats
path: root/buch/papers/laguerre/scripts
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--buch/papers/laguerre/scripts/estimates.py49
-rw-r--r--buch/papers/laguerre/scripts/gamma_approx.py116
-rw-r--r--buch/papers/laguerre/scripts/integrand.py42
-rw-r--r--buch/papers/laguerre/scripts/integrand_exp.py46
-rw-r--r--buch/papers/laguerre/scripts/laguerre_poly.py108
-rw-r--r--buch/papers/laguerre/scripts/rel_error_complex.py43
-rw-r--r--buch/papers/laguerre/scripts/rel_error_mirror.py38
-rw-r--r--buch/papers/laguerre/scripts/rel_error_range.py41
-rw-r--r--buch/papers/laguerre/scripts/rel_error_shifted.py40
-rw-r--r--buch/papers/laguerre/scripts/rel_error_simple.py41
-rw-r--r--buch/papers/laguerre/scripts/targets.py58
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")