From 766fd00f58d8aef648b447c908fc2c9cad83cbc5 Mon Sep 17 00:00:00 2001 From: Valentin Boettcher Date: Sat, 23 Jul 2022 18:05:58 +0200 Subject: [PATCH] narrow down roots in old gaussflow as well --- hopsflow/gaussflow.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/hopsflow/gaussflow.py b/hopsflow/gaussflow.py index 38d8d87..02b5d5a 100644 --- a/hopsflow/gaussflow.py +++ b/hopsflow/gaussflow.py @@ -6,6 +6,7 @@ import numpy.typing as npt from . import util from numpy.polynomial import Polynomial from .util import BCF, expand_t +import scipy @dataclass @@ -34,6 +35,15 @@ class SystemParams: W: np.ndarray = field(init=False) G: np.ndarray = field(init=False) + root_tol: float = 1e-7 + """ + The relative tolerance of the roots. + + The roots are being calculated with the numpy polynomial module + and then refined with newtons method. This parameter is passed as + ``rtol`` to :any:`scipy.optimize.newton`. + """ + def __post_init__(self): self.t_max = self.α_0.t_max @@ -115,6 +125,20 @@ def calculate_coefficients(sys: SystemParams) -> tuple[np.ndarray, np.ndarray]: You can try to alter the number of terms in the expansion of the BCF.""" ) + # improving the accuracy of the roots seems to be essential + p_prime = p.deriv() + p_prime_2 = p_prime.deriv() + for i, root in enumerate(master_roots): + res = scipy.optimize.newton( + p, + root, + p_prime, + rtol=sys.root_tol, + fprime2=p_prime_2, + maxiter=100000, # we can allow that, there aren't that many roots + ) + master_roots[i] = res + resiquals = f_0(master_roots) / p.deriv()(master_roots) return master_roots, resiquals