21#ifndef OPM_ROOTFINDERS_HEADER
22#define OPM_ROOTFINDERS_HEADER
24#include <opm/common/ErrorMacros.hpp>
25#include <opm/common/OpmLog/OpmLog.hpp>
38 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
40 std::ostringstream sstr;
41 sstr <<
"Error in parameters, zero not bracketed: [a, b] = ["
42 << x0 <<
", " << x1 <<
"] f(a) = " << f0 <<
" f(b) = " << f1;
43 OpmLog::debug(sstr.str());
44 OPM_THROW_NOLOG(std::runtime_error, sstr.str());
47 static double handleTooManyIterations(
const double x0,
const double x1,
const int maxiter)
49 OPM_THROW(std::runtime_error,
"Maximum number of iterations exceeded: " << maxiter <<
"\n"
50 <<
"Current interval is [" << std::min(x0, x1) <<
", "
51 << std::max(x0, x1) <<
"] abs(x0-x1) " << std::abs(x0-x1));
59 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
62 std::cerr <<
"Error in parameters, zero not bracketed: [a, b] = ["
63 << x0 <<
", " << x1 <<
"] f(a) = " << f0 <<
" f(b) = " << f1
65 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
67 static double handleTooManyIterations(
const double x0,
const double x1,
const int maxiter)
70 std::cerr <<
"Maximum number of iterations exceeded: " << maxiter
71 <<
", current interval is [" << std::min(x0, x1) <<
", "
72 << std::max(x0, x1) <<
"] abs(x0-x1) " << std::abs(x0-x1);
80 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
82 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
84 static double handleTooManyIterations(
const double x0,
const double x1,
const int )
92 template <
class ErrorPolicy = ThrowOnError>
103 template <
class Functor>
104 inline static double solve(
const Functor& f,
108 const double tolerance,
109 int& iterations_used)
112 const double macheps = numeric_limits<double>::epsilon();
113 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
118 const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
119 if (fabs(f0) < epsF) {
123 if (fabs(f1) < epsF) {
127 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
132 while (fabs(x1 - x0) >= eps) {
133 double xnew = regulaFalsiStep(x0, x1, f0, f1);
134 double fnew = f(xnew);
137 if (iterations_used > max_iter) {
138 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
140 if (fabs(fnew) < epsF) {
144 if ((fnew > 0.0) == (f0 > 0.0)) {
166 const double gamma = f1/(f1 + fnew);
172 return 0.5*(x0 + x1);
182 template <
class Functor>
183 inline static double solve(
const Functor& f,
184 const double initial_guess,
188 const double tolerance,
189 int& iterations_used)
192 const double macheps = numeric_limits<double>::epsilon();
193 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
195 double f_initial = f(initial_guess);
196 const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
197 if (fabs(f_initial) < epsF) {
198 return initial_guess;
202 double f0 = f_initial;
203 double f1 = f_initial;
204 if (x0 != initial_guess) {
206 if (fabs(f0) < epsF) {
210 if (x1 != initial_guess) {
212 if (fabs(f1) < epsF) {
216 if (f0*f_initial < 0.0) {
224 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
229 while (fabs(x1 - x0) >= eps) {
230 double xnew = regulaFalsiStep(x0, x1, f0, f1);
231 double fnew = f(xnew);
234 if (iterations_used > max_iter) {
235 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
237 if (fabs(fnew) < epsF) {
241 if ((fnew > 0.0) == (f0 > 0.0)) {
263 const double gamma = f1/(f1 + fnew);
269 return 0.5*(x0 + x1);
274 inline static double regulaFalsiStep(
const double a,
280 return (b*fa - a*fb)/(fa - fb);
288 template <
class ErrorPolicy = ThrowOnError>
292 inline static std::string name()
294 return "RegulaFalsiBisection";
303 template <
class Functor>
304 inline static double solve(
const Functor& f,
308 const double tolerance,
309 int& iterations_used)
312 const double sqrt_half = std::sqrt(0.5);
313 const double macheps = numeric_limits<double>::epsilon();
314 const double eps = tolerance + macheps * max(max(fabs(a), fabs(b)), 1.0);
319 const double epsF = tolerance + macheps * max(fabs(f0), 1.0);
320 if (fabs(f0) < epsF) {
324 if (fabs(f1) < epsF) {
328 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
333 double width = fabs(x1 - x0);
334 double contraction = 1.0;
335 while (width >= eps) {
340 const double xnew = (contraction < sqrt_half)
341 ? regulaFalsiStep(x0, x1, f0, f1)
342 : bisectionStep(x0, x1, f0, f1);
343 const double fnew = f(xnew);
345 if (iterations_used > max_iter) {
346 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
348 if (fabs(fnew) < epsF) {
352 if ((fnew > 0.0) == (f0 > 0.0)) {
374 const double gamma = f1 / (f1 + fnew);
381 const double widthnew = fabs(x1 - x0);
382 contraction = widthnew/width;
385 return 0.5 * (x0 + x1);
390 inline static double regulaFalsiStep(
const double a,
const double b,
const double fa,
const double fb)
392 assert(fa * fb < 0.0);
393 return (b * fa - a * fb) / (fa - fb);
395 inline static double bisectionStep(
const double a,
const double b,
const double fa,
const double fb)
397 static_cast<void>(fa);
398 static_cast<void>(fb);
399 assert(fa * fb < 0.0);
410 template <
class Functor>
417 const int max_iters = 100;
421 for (; i < max_iters; ++i) {
422 double x = x0 + cur_dx;
424 if (f0*f_new <= 0.0) {
427 cur_dx = -2.0*cur_dx;
429 if (i == max_iters) {
430 OPM_THROW(std::runtime_error,
"Could not bracket zero in " << max_iters <<
"iterations.");
434 b = i < 2 ? x0 : x0 + 0.25*cur_dx;
436 a = i < 2 ? x0 : x0 + 0.25*cur_dx;
Definition: RootFinders.hpp:290
static double solve(const Functor &f, const double a, const double b, const int max_iter, const double tolerance, int &iterations_used)
Implements a modified regula falsi method as described in "Improved algorithms of Illinois-type for t...
Definition: RootFinders.hpp:304
Definition: RootFinders.hpp:94
static double solve(const Functor &f, const double initial_guess, const double a, const double b, const int max_iter, const double tolerance, int &iterations_used)
Implements a modified regula falsi method as described in "Improved algorithms of Illinois-type for t...
Definition: RootFinders.hpp:183
static double solve(const Functor &f, const double a, const double b, const int max_iter, const double tolerance, int &iterations_used)
Implements a modified regula falsi method as described in "Improved algorithms of Illinois-type for t...
Definition: RootFinders.hpp:104
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:29
void bracketZero(const Functor &f, const double x0, const double dx, double &a, double &b)
Attempts to find an interval bracketing a zero by successive enlargement of search interval.
Definition: RootFinders.hpp:411
Definition: RootFinders.hpp:79
Definition: RootFinders.hpp:37
Definition: RootFinders.hpp:58