21#ifndef OPM_ROOTFINDERS_HEADER
22#define OPM_ROOTFINDERS_HEADER
24#include <opm/common/ErrorMacros.hpp>
35 static double handleBracketingFailure(
const double x0,
const double x1,
36 const double f0,
const double f1);
38 static double handleTooManyIterations(
const double x0,
39 const double x1,
const int maxiter);
45 static double handleBracketingFailure(
const double x0,
const double x1,
46 const double f0,
const double f1);
47 static double handleTooManyIterations(
const double x0,
48 const double x1,
const int maxiter);
54 static double handleBracketingFailure(
const double x0,
const double x1,
55 const double f0,
const double f1)
57 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
59 static double handleTooManyIterations(
const double x0,
60 const double x1,
const int )
68 template <
class ErrorPolicy = ThrowOnError>
79 template <
class Functor>
80 inline static double solve(
const Functor& f,
84 const double tolerance,
88 const double macheps = numeric_limits<double>::epsilon();
89 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
94 const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
95 if (fabs(f0) < epsF) {
99 if (fabs(f1) < epsF) {
103 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
108 while (fabs(x1 - x0) >= eps) {
109 double xnew = regulaFalsiStep(x0, x1, f0, f1);
110 double fnew = f(xnew);
113 if (iterations_used > max_iter) {
114 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
116 if (fabs(fnew) < epsF) {
120 if ((fnew > 0.0) == (f0 > 0.0)) {
142 const double gamma = f1/(f1 + fnew);
148 return 0.5*(x0 + x1);
158 template <
class Functor>
159 inline static double solve(
const Functor& f,
160 const double initial_guess,
164 const double tolerance,
165 int& iterations_used)
168 const double macheps = numeric_limits<double>::epsilon();
169 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
171 double f_initial = f(initial_guess);
172 const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
173 if (fabs(f_initial) < epsF) {
174 return initial_guess;
178 double f0 = f_initial;
179 double f1 = f_initial;
180 if (x0 != initial_guess) {
182 if (fabs(f0) < epsF) {
186 if (x1 != initial_guess) {
188 if (fabs(f1) < epsF) {
192 if (f0*f_initial < 0.0) {
200 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
205 while (fabs(x1 - x0) >= eps) {
206 double xnew = regulaFalsiStep(x0, x1, f0, f1);
207 double fnew = f(xnew);
210 if (iterations_used > max_iter) {
211 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
213 if (fabs(fnew) < epsF) {
217 if ((fnew > 0.0) == (f0 > 0.0)) {
239 const double gamma = f1/(f1 + fnew);
245 return 0.5*(x0 + x1);
250 inline static double regulaFalsiStep(
const double a,
256 return (b*fa - a*fb)/(fa - fb);
264 template <
class ErrorPolicy = ThrowOnError>
268 inline static std::string name()
270 return "RegulaFalsiBisection";
279 template <
class Functor>
280 inline static double solve(
const Functor& f,
284 const double tolerance,
285 int& iterations_used)
288 const double sqrt_half = std::sqrt(0.5);
289 const double macheps = numeric_limits<double>::epsilon();
290 const double eps = tolerance + macheps * max(max(fabs(a), fabs(b)), 1.0);
295 const double epsF = tolerance + macheps * max(fabs(f0), 1.0);
296 if (fabs(f0) < epsF) {
300 if (fabs(f1) < epsF) {
304 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
309 double width = fabs(x1 - x0);
310 double contraction = 1.0;
311 while (width >= eps) {
316 const double xnew = (contraction < sqrt_half)
317 ? regulaFalsiStep(x0, x1, f0, f1)
318 : bisectionStep(x0, x1, f0, f1);
319 const double fnew = f(xnew);
321 if (iterations_used > max_iter) {
322 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
324 if (fabs(fnew) < epsF) {
328 if ((fnew > 0.0) == (f0 > 0.0)) {
350 const double gamma = f1 / (f1 + fnew);
357 const double widthnew = fabs(x1 - x0);
358 contraction = widthnew/width;
361 return 0.5 * (x0 + x1);
366 inline static double regulaFalsiStep(
const double a,
const double b,
const double fa,
const double fb)
368 assert(fa * fb < 0.0);
369 return (b * fa - a * fb) / (fa - fb);
371 inline static double bisectionStep(
const double a,
const double b,
const double fa,
const double fb)
373 static_cast<void>(fa);
374 static_cast<void>(fb);
375 assert(fa * fb < 0.0);
386 template <
class Functor>
393 const int max_iters = 100;
397 for (; i < max_iters; ++i) {
398 double x = x0 + cur_dx;
400 if (f0*f_new <= 0.0) {
403 cur_dx = -2.0*cur_dx;
405 if (i == max_iters) {
406 OPM_THROW(std::runtime_error,
407 "Could not bracket zero in " +
408 std::to_string(max_iters) +
" iterations.");
412 b = i < 2 ? x0 : x0 + 0.25*cur_dx;
414 a = i < 2 ? x0 : x0 + 0.25*cur_dx;
Definition RootFinders.hpp:266
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:280
Definition RootFinders.hpp:70
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:159
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:80
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
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:387
Definition RootFinders.hpp:53
Definition RootFinders.hpp:34
Definition RootFinders.hpp:44