Source code for triqs.utility.bound_and_bisect

# Copyright (c) 2014-2016 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
# Copyright (c) 2014-2016 Centre national de la recherche scientifique (CNRS)
# Copyright (c) 2020 Simons Foundation
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You may obtain a copy of the License at
#     https:#www.gnu.org/licenses/gpl-3.0.txt
#
# Authors: Priyanka Seth, Nils Wentzell

import scipy.optimize

[docs]def determine_bounds(F, x_0, dx, maxiter): """ Determine the bounds [a,b] for the zero of function F, using x_0 as a guess and incrementing by dx each time.""" # Determine if function is increasing or decreasing if F(x_0 + dx) - F(x_0) < 0: dx *= -1.0 # Determine lower bound x = x_0; i = 0 while not F(x) < 0.0: x -= dx i += 1 if i > maxiter: raise ValueError("determine_bounds: cannot find lower bound.") a = x # Determine upper bound x = x_0; i = 0 while not F(x) > 0.0: x += dx i += 1 if i > maxiter: raise ValueError("determine_bounds: cannot find upper bound.") b = x return a, b
[docs]def bound_and_bisect(f, x_0, y=0.0, dx=1.0, xtol=1e-3, x_name='x', y_name='y', maxiter=1000, verbosity=1): """ Solve F(x) = y where F is a monotonic function. Arguments : - F: function (real valued) to be solved - x_0: Initial guess for for x. - y: If given, F(x) = y will be solved. - dx: Increment to use in bound search. - x_name, y_name: Names of the variables x and y for the report. - maxiter: Maximum number of iterations for both bound searching and bisection. Returns : - A tuple (x,f(x)). """ F = lambda x: f(x) - y a,b = determine_bounds(F,x_0,dx,maxiter) if verbosity > 0: print('Bounds are determined: %s <= %s <= %s'%(a,x_name,b)) x = scipy.optimize.bisect(F, a, b, xtol=xtol, maxiter=maxiter) fx = f(x) if verbosity > 0: print('%s(%s) = %s solved:'%(y_name,x_name,y)) print('%s = %s => %s = %s'%(x_name, x, y_name, fx)) return x,fx