Searching for Roots

std::vector<Root> DCProgs::find_roots(DeterminantEq const &_det, t_real _xtol = 1e-8, t_real _rtol = 1e-8, t_uint _itermax = 100, t_real _lowerbound = quiet_nan, t_real _upperbound = quiet_nan)

Finds root using brentq and find_root_intervals.

Tries and computes the roots of an input determinant equation. This is a three fold process:

  1. Find an interval that contains all roots/eigenvalues
  1. Split it into smaller intervals until each contains only one root
  1. Optimize over this interval to find an accurate position for the rootThe last step is carried out by brentq().

Parameters
  • _det -

    The determinant equation for which to compute the roots.

  • _xtol -

    Tolerance for interval size

  • _rtol -

    Tolerance for interval size. The convergence criteria is an affine function of the root: \(x_{\mathrm{tol}} + r_{\mathrm{tol}} x_{\mathrm{current}} = \frac{|x_a - x_b|}{2}\).

  • _itermax -

    maximum number of iterations for any of the three steps.

  • _lowerbound -

    Lower bound of the interval bracketing all roots. If None, the lower bound is obtained from find_lower_bound_for_roots().

  • _upperbound -

    Upper bound of the interval bracketing all roots. If None, the upper bound is obtained from find_upper_bound_for_roots().

Bracketing all Roots

t_real DCProgs::find_lower_bound_for_roots(DeterminantEq const &_det, t_real _start = 0e0, t_real _alpha = 5e0, t_uint _itermax = 100)

Figures out an lower bound for root finding.

Proceeds by computing the eigenvalues, then setting lower bound to somewhat lower than the lowest eigenvalue. It then checks that the eigenvalues of the matrix computed at that value, and so on and so forth. The algorithm stops when the lowest eigenvalue is higher than the current bound.

Parameters
  • _det -

    The determinant equation

  • _start -

    Value where to start looking for lower bound.

  • _alpha -

    factor by which to set new lower bound: \(s_{n+1} = min(\epsilon_i) + \alpha (s_N - min(\epsilon_i))\).

  • _itermax -

    Maximum number of iterations.

t_real DCProgs::find_upper_bound_for_roots(DeterminantEq const &_det, t_real _start = 0e0, t_real _alpha = 5e0, t_uint _itermax = 100)

Figures out an upper bound for root finding.

Parameters
  • _det -

    The determinant equation

  • _start -

    Value where to start looking for lower bound.

  • _alpha -

    factor by which to set new lower bound: \(s_{n+1} = min(\epsilon_i) + \alpha (s_N - min(\epsilon_i))\).

  • _itermax -

    Maximum number of iterations.

Finding intervals for each root

std::vector<RootInterval> DCProgs::find_root_intervals(DeterminantEq const &_det, t_real _mins = quiet_nan, t_real _maxs = quiet_nan, t_real _tolerance = 1e-8)

Figures out interval where roots can be found.

Parameters
  • _det -

    The determinant equation

  • _mins -

    A valid lower bound, or DCProgs::quiet_nan. In the latter case, the lower bound is determined using DCProgs::find_lower_bound_for_roots().

  • _maxs -

    A valid upper bound, or DCProgs::quiet_nan. In the latter case, the upper bound is determined using DCProgs::find_upper_bound_for_roots().

  • _tolerance -

    Minimum size of intervals. Below that, roots are expected to be multiples.

std::vector<RootInterval> DCProgs::find_root_intervals_brute_force(DeterminantEq const &_det, t_real _resolution = 1e-1, t_real _mins = quiet_nan, t_real _maxs = quiet_nan, t_real _tolerance = 1e-1)

Finds roots via brute force search.

Computes all values between mins and maxs, for a given resolution. If determinant changes sign between two values, or if it comes to within tolerance of zero, then computes eigenvalues of H to determine possible multiplicity.

Parameters
  • _det -

    The determinant equation

  • _resolution -

    resolution at which computes values in interval.

  • _mins -

    A valid lower bound, or DCProgs::quiet_nan. In the latter case, the lower bound is determined using DCProgs::find_lower_bound_for_roots().

  • _maxs -

    A valid upper bound, or DCProgs::quiet_nan. In the latter case, the upper bound is determined using DCProgs::find_upper_bound_for_roots().

  • _tolerance -

    Tolerance below which the value of the determinant is considered “close to zero”.

Root optimization

std::tuple<t_real, t_uint, t_uint> DCProgs::brentq(std::function<t_real(t_real)> const &_function, t_real _xstart, t_real _xend, t_real _xtol = 1e-8, t_real _rtol = 1e-8, t_uint _itermax = 100, )

Computes root of a function in a given interval.

Scavenged from Scipy. Actual code (.cc file) is under BSD.

Return
A tuple (x, iterations, function calls)
Parameters
  • _function -

    A call-back to the actual function.

  • _xstart -

    Beginning of the interval

  • _xend -

    End of the interval

  • _xtol -

    Tolerance for interval size

  • _rtol -

    Tolerance for interval size. The convergence criteria is an affine function of the root: \(x_{\mathrm{tol}} + r_{\mathrm{tol}} x_{\mathrm{current}} = \frac{|x_a - x_b|}{2}\).

  • _itermax -

    maximum number of iterations.