diff options
Diffstat (limited to 'buch/papers/laguerre/scripts/gamma_approx.ipynb')
| -rw-r--r-- | buch/papers/laguerre/scripts/gamma_approx.ipynb | 80 |
1 files changed, 68 insertions, 12 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": { |