From 3cb2fa354f814fa98474610dac744281285dafc6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Fri, 15 Jul 2022 11:40:55 +0200
Subject: First version of section 'Gauss Quadratur', fix to gamma_approx.py
 when z=0

---
 buch/papers/laguerre/scripts/gamma_approx.py    | 18 +++++++++++++-----
 buch/papers/laguerre/scripts/rel_error_range.py |  2 +-
 2 files changed, 14 insertions(+), 6 deletions(-)

(limited to 'buch/papers/laguerre/scripts')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 208f770..9f9dae7 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -1,7 +1,5 @@
 from pathlib import Path
 
-import matplotlib as mpl
-import matplotlib.pyplot as plt
 import numpy as np
 import scipy.special
 
@@ -58,6 +56,8 @@ def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
 
 
 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)
 
@@ -73,6 +73,8 @@ def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
 
 
 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)
@@ -81,6 +83,8 @@ def laguerre_gamma_simple(z, x=None, w=None, n=8):
 
 
 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:
@@ -90,8 +94,8 @@ def laguerre_gamma_mirror(z, x=None, w=None, n=8):
     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)
+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":
@@ -104,4 +108,8 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
 
 
 def calc_rel_error(x, y):
-    return (y - x) / x
+    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/rel_error_range.py b/buch/papers/laguerre/scripts/rel_error_range.py
index 7d017a7..7c74d76 100644
--- a/buch/papers/laguerre/scripts/rel_error_range.py
+++ b/buch/papers/laguerre/scripts/rel_error_range.py
@@ -5,7 +5,7 @@ if __name__ == "__main__":
 
     import gamma_approx as ga
 
-    N = 1000
+    N = 1001
     xmin = -5
     xmax = 5
     ns = np.arange(2, 12, 2)
-- 
cgit v1.2.1