My Project
RootFinders.hpp
1 /*
2  Copyright 2010, 2019 SINTEF Digital
3  Copyright 2010, 2019 Equinor ASA
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_ROOTFINDERS_HEADER
22 #define OPM_ROOTFINDERS_HEADER
23 
24 #include <opm/common/ErrorMacros.hpp>
25 #include <opm/common/OpmLog/OpmLog.hpp>
26 
27 #include <string>
28 #include <algorithm>
29 #include <limits>
30 #include <cmath>
31 #include <iostream>
32 
33 namespace Opm
34 {
35 
36  struct ThrowOnError
37  {
38  static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
39  {
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());
45  return -1e100; // Never reached.
46  }
47  static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
48  {
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));
52  return -1e100; // Never reached.
53  }
54  };
55 
56 
58  {
59  static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
60  {
61  OPM_REPORT;
62  std::cerr << "Error in parameters, zero not bracketed: [a, b] = ["
63  << x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1
64  << "";
65  return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
66  }
67  static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
68  {
69  OPM_REPORT;
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);
73  return 0.5*(x0 + x1);
74  }
75  };
76 
77 
79  {
80  static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
81  {
82  return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
83  }
84  static double handleTooManyIterations(const double x0, const double x1, const int /*maxiter*/)
85  {
86  return 0.5*(x0 + x1);
87  }
88  };
89 
90 
91 
92  template <class ErrorPolicy = ThrowOnError>
94  {
95  public:
96 
97 
103  template <class Functor>
104  inline static double solve(const Functor& f,
105  const double a,
106  const double b,
107  const int max_iter,
108  const double tolerance,
109  int& iterations_used)
110  {
111  using namespace std;
112  const double macheps = numeric_limits<double>::epsilon();
113  const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
114 
115  double x0 = a;
116  double x1 = b;
117  double f0 = f(x0);
118  const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
119  if (fabs(f0) < epsF) {
120  return x0;
121  }
122  double f1 = f(x1);
123  if (fabs(f1) < epsF) {
124  return x1;
125  }
126  if (f0*f1 > 0.0) {
127  return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
128  }
129  iterations_used = 0;
130  // In every iteraton, x1 is the last point computed,
131  // and x0 is the last point computed that makes it a bracket.
132  while (fabs(x1 - x0) >= eps) {
133  double xnew = regulaFalsiStep(x0, x1, f0, f1);
134  double fnew = f(xnew);
135 // cout << "xnew = " << xnew << " fnew = " << fnew << endl;
136  ++iterations_used;
137  if (iterations_used > max_iter) {
138  return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
139  }
140  if (fabs(fnew) < epsF) {
141  return xnew;
142  }
143  // Now we must check which point we must replace.
144  if ((fnew > 0.0) == (f0 > 0.0)) {
145  // We must replace x0.
146  x0 = x1;
147  f0 = f1;
148  } else {
149  // We must replace x1, this is the case where
150  // the modification to regula falsi kicks in,
151  // by modifying f0.
152  // 1. The classic Illinois method
153 // const double gamma = 0.5;
154  // @afr: The next two methods do not work??!!?
155  // 2. The method called 'Method 3' in the paper.
156 // const double phi0 = f1/f0;
157 // const double phi1 = fnew/f1;
158 // const double gamma = 1.0 - phi1/(1.0 - phi0);
159  // 3. The method called 'Method 4' in the paper.
160 // const double phi0 = f1/f0;
161 // const double phi1 = fnew/f1;
162 // const double gamma = 1.0 - phi0 - phi1;
163 // cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
164 // " gamma = " << gamma << endl;
165  // 4. The 'Pegasus' method
166  const double gamma = f1/(f1 + fnew);
167  f0 *= gamma;
168  }
169  x1 = xnew;
170  f1 = fnew;
171  }
172  return 0.5*(x0 + x1);
173  }
174 
175 
182  template <class Functor>
183  inline static double solve(const Functor& f,
184  const double initial_guess,
185  const double a,
186  const double b,
187  const int max_iter,
188  const double tolerance,
189  int& iterations_used)
190  {
191  using namespace std;
192  const double macheps = numeric_limits<double>::epsilon();
193  const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
194 
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;
199  }
200  double x0 = a;
201  double x1 = b;
202  double f0 = f_initial;
203  double f1 = f_initial;
204  if (x0 != initial_guess) {
205  f0 = f(x0);
206  if (fabs(f0) < epsF) {
207  return x0;
208  }
209  }
210  if (x1 != initial_guess) {
211  f1 = f(x1);
212  if (fabs(f1) < epsF) {
213  return x1;
214  }
215  }
216  if (f0*f_initial < 0.0) {
217  x1 = initial_guess;
218  f1 = f_initial;
219  } else {
220  x0 = initial_guess;
221  f0 = f_initial;
222  }
223  if (f0*f1 > 0.0) {
224  return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
225  }
226  iterations_used = 0;
227  // In every iteraton, x1 is the last point computed,
228  // and x0 is the last point computed that makes it a bracket.
229  while (fabs(x1 - x0) >= eps) {
230  double xnew = regulaFalsiStep(x0, x1, f0, f1);
231  double fnew = f(xnew);
232 // cout << "xnew = " << xnew << " fnew = " << fnew << endl;
233  ++iterations_used;
234  if (iterations_used > max_iter) {
235  return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
236  }
237  if (fabs(fnew) < epsF) {
238  return xnew;
239  }
240  // Now we must check which point we must replace.
241  if ((fnew > 0.0) == (f0 > 0.0)) {
242  // We must replace x0.
243  x0 = x1;
244  f0 = f1;
245  } else {
246  // We must replace x1, this is the case where
247  // the modification to regula falsi kicks in,
248  // by modifying f0.
249  // 1. The classic Illinois method
250 // const double gamma = 0.5;
251  // @afr: The next two methods do not work??!!?
252  // 2. The method called 'Method 3' in the paper.
253 // const double phi0 = f1/f0;
254 // const double phi1 = fnew/f1;
255 // const double gamma = 1.0 - phi1/(1.0 - phi0);
256  // 3. The method called 'Method 4' in the paper.
257 // const double phi0 = f1/f0;
258 // const double phi1 = fnew/f1;
259 // const double gamma = 1.0 - phi0 - phi1;
260 // cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
261 // " gamma = " << gamma << endl;
262  // 4. The 'Pegasus' method
263  const double gamma = f1/(f1 + fnew);
264  f0 *= gamma;
265  }
266  x1 = xnew;
267  f1 = fnew;
268  }
269  return 0.5*(x0 + x1);
270  }
271 
272 
273  private:
274  inline static double regulaFalsiStep(const double a,
275  const double b,
276  const double fa,
277  const double fb)
278  {
279  assert(fa*fb < 0.0);
280  return (b*fa - a*fb)/(fa - fb);
281  }
282 
283 
284  };
285 
286 
287 
288  template <class ErrorPolicy = ThrowOnError>
290  {
291  public:
292  inline static std::string name()
293  {
294  return "RegulaFalsiBisection";
295  }
296 
303  template <class Functor>
304  inline static double solve(const Functor& f,
305  const double a,
306  const double b,
307  const int max_iter,
308  const double tolerance,
309  int& iterations_used)
310  {
311  using namespace std;
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);
315 
316  double x0 = a;
317  double x1 = b;
318  double f0 = f(x0);
319  const double epsF = tolerance + macheps * max(fabs(f0), 1.0);
320  if (fabs(f0) < epsF) {
321  return x0;
322  }
323  double f1 = f(x1);
324  if (fabs(f1) < epsF) {
325  return x1;
326  }
327  if (f0 * f1 > 0.0) {
328  return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
329  }
330  iterations_used = 0;
331  // In every iteraton, x1 is the last point computed,
332  // and x0 is the last point computed that makes it a bracket.
333  double width = fabs(x1 - x0);
334  double contraction = 1.0;
335  while (width >= eps) {
336  // If we are contracting sufficiently to at least halve
337  // the interval width in two iterations we use regula
338  // falsi. Otherwise, we take a bisection step to avoid
339  // slow convergence.
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);
344  ++iterations_used;
345  if (iterations_used > max_iter) {
346  return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
347  }
348  if (fabs(fnew) < epsF) {
349  return xnew;
350  }
351  // Now we must check which point we must replace.
352  if ((fnew > 0.0) == (f0 > 0.0)) {
353  // We must replace x0.
354  x0 = x1;
355  f0 = f1;
356  } else {
357  // We must replace x1, this is the case where
358  // the modification to regula falsi kicks in,
359  // by modifying f0.
360  // 1. The classic Illinois method
361  // const double gamma = 0.5;
362  // @afr: The next two methods do not work??!!?
363  // 2. The method called 'Method 3' in the paper.
364  // const double phi0 = f1/f0;
365  // const double phi1 = fnew/f1;
366  // const double gamma = 1.0 - phi1/(1.0 - phi0);
367  // 3. The method called 'Method 4' in the paper.
368  // const double phi0 = f1/f0;
369  // const double phi1 = fnew/f1;
370  // const double gamma = 1.0 - phi0 - phi1;
371  // cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
372  // " gamma = " << gamma << endl;
373  // 4. The 'Pegasus' method
374  const double gamma = f1 / (f1 + fnew);
375  // 5. Straight unmodified Regula Falsi
376  // const double gamma = 1.0;
377  f0 *= gamma;
378  }
379  x1 = xnew;
380  f1 = fnew;
381  const double widthnew = fabs(x1 - x0);
382  contraction = widthnew/width;
383  width = widthnew;
384  }
385  return 0.5 * (x0 + x1);
386  }
387 
388 
389  private:
390  inline static double regulaFalsiStep(const double a, const double b, const double fa, const double fb)
391  {
392  assert(fa * fb < 0.0);
393  return (b * fa - a * fb) / (fa - fb);
394  }
395  inline static double bisectionStep(const double a, const double b, const double fa, const double fb)
396  {
397  static_cast<void>(fa);
398  static_cast<void>(fb);
399  assert(fa * fb < 0.0);
400  return 0.5*(a + b);
401  }
402  };
403 
404 
405 
406 
407 
410  template <class Functor>
411  inline void bracketZero(const Functor& f,
412  const double x0,
413  const double dx,
414  double& a,
415  double& b)
416  {
417  const int max_iters = 100;
418  double f0 = f(x0);
419  double cur_dx = dx;
420  int i = 0;
421  for (; i < max_iters; ++i) {
422  double x = x0 + cur_dx;
423  double f_new = f(x);
424  if (f0*f_new <= 0.0) {
425  break;
426  }
427  cur_dx = -2.0*cur_dx;
428  }
429  if (i == max_iters) {
430  OPM_THROW(std::runtime_error, "Could not bracket zero in " << max_iters << "iterations.");
431  }
432  if (cur_dx < 0.0) {
433  a = x0 + cur_dx;
434  b = i < 2 ? x0 : x0 + 0.25*cur_dx;
435  } else {
436  a = i < 2 ? x0 : x0 + 0.25*cur_dx;
437  b = x0 + cur_dx;
438  }
439  }
440 
441 
442 
443 
444 
445 } // namespace Opm
446 
447 
448 
449 
450 #endif // OPM_ROOTFINDERS_HEADER
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