Go to the documentation of this file.
21#ifndef OPM_ROOTFINDERS_HEADER
22#define OPM_ROOTFINDERS_HEADER
40 std::ostringstream sstr;
41 sstr << "Error in parameters, zero not bracketed: [a, b] = ["
42 << x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1;
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));
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;
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);
82 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
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>
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;
#define OPM_THROW(Exception, message) Definition: ErrorMacros.hpp:52
#define OPM_THROW_NOLOG(Exception, message) Definition: ErrorMacros.hpp:63
#define OPM_REPORT Definition: ErrorMacros.hpp:39
const cJSON *const b Definition: cJSON.h:251
const char *const string Definition: cJSON.h:170
static void debug(const std::string &message)
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) Definition: RootFinders.hpp:304
static std::string name() Definition: RootFinders.hpp:292
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) 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) Definition: RootFinders.hpp:104
void bracketZero(const Functor &f, const double x0, const double dx, double &a, double &b) Definition: RootFinders.hpp:411
T min(const T v0, const T v1) Definition: exprtk.hpp:1400
T max(const T v0, const T v1) Definition: exprtk.hpp:1407
x y t t *t x y t t t x y t t t x *y t *t t x *y t *t t x y t t t x y t t t x(y+z)
Definition: RootFinders.hpp:79
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) Definition: RootFinders.hpp:80
static double handleTooManyIterations(const double x0, const double x1, const int) Definition: RootFinders.hpp:84
Definition: RootFinders.hpp:37
static double handleTooManyIterations(const double x0, const double x1, const int maxiter) Definition: RootFinders.hpp:47
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) Definition: RootFinders.hpp:38
Definition: RootFinders.hpp:58
static double handleTooManyIterations(const double x0, const double x1, const int maxiter) Definition: RootFinders.hpp:67
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) Definition: RootFinders.hpp:59
|