36 #ifndef OPM_ROOTFINDERS_HEADER 37 #define OPM_ROOTFINDERS_HEADER 39 #include <opm/common/ErrorMacros.hpp> 51 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
53 OPM_THROW(std::runtime_error,
"Error in parameters, zero not bracketed: [a, b] = [" 54 << x0 <<
", " << x1 <<
"] f(a) = " << f0 <<
" f(b) = " << f1);
57 static double handleTooManyIterations(
const double x0,
const double x1,
const int maxiter)
59 OPM_THROW(std::runtime_error,
"Maximum number of iterations exceeded: " << maxiter <<
"\n" 60 <<
"Current interval is [" << std::min(x0, x1) <<
", " 61 << std::max(x0, x1) <<
"]");
69 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
72 std::cerr <<
"Error in parameters, zero not bracketed: [a, b] = [" 73 << x0 <<
", " << x1 <<
"] f(a) = " << f0 <<
" f(b) = " << f1
75 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
77 static double handleTooManyIterations(
const double x0,
const double x1,
const int maxiter)
80 std::cerr <<
"Maximum number of iterations exceeded: " << maxiter
81 <<
", current interval is [" << std::min(x0, x1) <<
", " 82 << std::max(x0, x1) <<
"]";
90 static double handleBracketingFailure(
const double x0,
const double x1,
const double f0,
const double f1)
92 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
94 static double handleTooManyIterations(
const double x0,
const double x1,
const int )
102 template <
class ErrorPolicy = ThrowOnError>
113 template <
class Functor>
114 inline static double solve(
const Functor& f,
118 const double tolerance,
119 int& iterations_used)
122 const double macheps = numeric_limits<double>::epsilon();
123 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
128 const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
129 if (fabs(f0) < epsF) {
133 if (fabs(f1) < epsF) {
137 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
142 while (fabs(x1 - x0) >= 1e-9*eps) {
143 double xnew = regulaFalsiStep(x0, x1, f0, f1);
144 double fnew = f(xnew);
147 if (iterations_used > max_iter) {
148 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
150 if (fabs(fnew) < epsF) {
154 if ((fnew > 0.0) == (f0 > 0.0)) {
176 const double gamma = f1/(f1 + fnew);
182 return 0.5*(x0 + x1);
192 template <
class Functor>
193 inline static double solve(
const Functor& f,
194 const double initial_guess,
198 const double tolerance,
199 int& iterations_used)
202 const double macheps = numeric_limits<double>::epsilon();
203 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
205 double f_initial = f(initial_guess);
206 const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
207 if (fabs(f_initial) < epsF) {
208 return initial_guess;
212 double f0 = f_initial;
213 double f1 = f_initial;
214 if (x0 != initial_guess) {
216 if (fabs(f0) < epsF) {
220 if (x1 != initial_guess) {
222 if (fabs(f1) < epsF) {
226 if (f0*f_initial < 0.0) {
234 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
239 while (fabs(x1 - x0) >= 1e-9*eps) {
240 double xnew = regulaFalsiStep(x0, x1, f0, f1);
241 double fnew = f(xnew);
244 if (iterations_used > max_iter) {
245 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
247 if (fabs(fnew) < epsF) {
251 if ((fnew > 0.0) == (f0 > 0.0)) {
273 const double gamma = f1/(f1 + fnew);
279 return 0.5*(x0 + x1);
284 inline static double regulaFalsiStep(
const double a,
290 return (b*fa - a*fb)/(fa - fb);
300 template <
class Functor>
301 inline void bracketZero(
const Functor& f,
307 const int max_iters = 100;
311 for (; i < max_iters; ++i) {
312 double x = x0 + cur_dx;
314 if (f0*f_new <= 0.0) {
317 cur_dx = -2.0*cur_dx;
319 if (i == max_iters) {
320 OPM_THROW(std::runtime_error,
"Could not bracket zero in " << max_iters <<
"iterations.");
324 b = i < 2 ? x0 : x0 + 0.25*cur_dx;
326 a = i < 2 ? x0 : x0 + 0.25*cur_dx;
340 #endif // OPM_ROOTFINDERS_HEADER Definition: RootFinders.hpp:103
Definition: RootFinders.hpp:49
Definition: AnisotropicEikonal.cpp:446
Definition: RootFinders.hpp:88
Definition: RootFinders.hpp:67
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:114
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:193