diff options
author | Patrik Müller <patrik.mueller@ost.ch> | 2022-05-31 16:31:25 +0200 |
---|---|---|
committer | Patrik Müller <patrik.mueller@ost.ch> | 2022-05-31 16:31:25 +0200 |
commit | 6149839224755c21225d2decddeae12207c2cbab (patch) | |
tree | 3b7a2f02ad8e388c68eb6b5f3ea4144a50c5ad56 /buch/papers/laguerre/scripts/gamma_approx.py | |
parent | Merge branch 'AndreasFMueller:master' into master (diff) | |
download | SeminarSpezielleFunktionen-6149839224755c21225d2decddeae12207c2cbab.tar.gz SeminarSpezielleFunktionen-6149839224755c21225d2decddeae12207c2cbab.zip |
Add rule of thumb, analyse integrand, correct mistake in integration SLP<->LP
Diffstat (limited to 'buch/papers/laguerre/scripts/gamma_approx.py')
-rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.py | 197 |
1 files changed, 197 insertions, 0 deletions
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() |