My Project
Loading...
Searching...
No Matches
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
26#include <cmath>
27#include <limits>
28#include <string>
29
30namespace Opm
31{
32
34 {
35 static double handleBracketingFailure(const double x0, const double x1,
36 const double f0, const double f1);
37
38 static double handleTooManyIterations(const double x0,
39 const double x1, const int maxiter);
40 };
41
42
44 {
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);
49 };
50
51
53 {
54 static double handleBracketingFailure(const double x0, const double x1,
55 const double f0, const double f1)
56 {
57 return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
58 }
59 static double handleTooManyIterations(const double x0,
60 const double x1, const int /*maxiter*/)
61 {
62 return 0.5*(x0 + x1);
63 }
64 };
65
66
67
68 template <class ErrorPolicy = ThrowOnError>
70 {
71 public:
72
73
79 template <class Functor>
80 inline static double solve(const Functor& f,
81 const double a,
82 const double b,
83 const int max_iter,
84 const double tolerance,
85 int& iterations_used)
86 {
87 using namespace std;
88 const double macheps = numeric_limits<double>::epsilon();
89 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
90
91 double x0 = a;
92 double x1 = b;
93 double f0 = f(x0);
94 const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
95 if (fabs(f0) < epsF) {
96 return x0;
97 }
98 double f1 = f(x1);
99 if (fabs(f1) < epsF) {
100 return x1;
101 }
102 if (f0*f1 > 0.0) {
103 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
104 }
105 iterations_used = 0;
106 // In every iteraton, x1 is the last point computed,
107 // and x0 is the last point computed that makes it a bracket.
108 while (fabs(x1 - x0) >= eps) {
109 double xnew = regulaFalsiStep(x0, x1, f0, f1);
110 double fnew = f(xnew);
111// cout << "xnew = " << xnew << " fnew = " << fnew << endl;
112 ++iterations_used;
113 if (iterations_used > max_iter) {
114 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
115 }
116 if (fabs(fnew) < epsF) {
117 return xnew;
118 }
119 // Now we must check which point we must replace.
120 if ((fnew > 0.0) == (f0 > 0.0)) {
121 // We must replace x0.
122 x0 = x1;
123 f0 = f1;
124 } else {
125 // We must replace x1, this is the case where
126 // the modification to regula falsi kicks in,
127 // by modifying f0.
128 // 1. The classic Illinois method
129// const double gamma = 0.5;
130 // @afr: The next two methods do not work??!!?
131 // 2. The method called 'Method 3' in the paper.
132// const double phi0 = f1/f0;
133// const double phi1 = fnew/f1;
134// const double gamma = 1.0 - phi1/(1.0 - phi0);
135 // 3. The method called 'Method 4' in the paper.
136// const double phi0 = f1/f0;
137// const double phi1 = fnew/f1;
138// const double gamma = 1.0 - phi0 - phi1;
139// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
140// " gamma = " << gamma << endl;
141 // 4. The 'Pegasus' method
142 const double gamma = f1/(f1 + fnew);
143 f0 *= gamma;
144 }
145 x1 = xnew;
146 f1 = fnew;
147 }
148 return 0.5*(x0 + x1);
149 }
150
151
158 template <class Functor>
159 inline static double solve(const Functor& f,
160 const double initial_guess,
161 const double a,
162 const double b,
163 const int max_iter,
164 const double tolerance,
165 int& iterations_used)
166 {
167 using namespace std;
168 const double macheps = numeric_limits<double>::epsilon();
169 const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
170
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;
175 }
176 double x0 = a;
177 double x1 = b;
178 double f0 = f_initial;
179 double f1 = f_initial;
180 if (x0 != initial_guess) {
181 f0 = f(x0);
182 if (fabs(f0) < epsF) {
183 return x0;
184 }
185 }
186 if (x1 != initial_guess) {
187 f1 = f(x1);
188 if (fabs(f1) < epsF) {
189 return x1;
190 }
191 }
192 if (f0*f_initial < 0.0) {
193 x1 = initial_guess;
194 f1 = f_initial;
195 } else {
196 x0 = initial_guess;
197 f0 = f_initial;
198 }
199 if (f0*f1 > 0.0) {
200 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
201 }
202 iterations_used = 0;
203 // In every iteraton, x1 is the last point computed,
204 // and x0 is the last point computed that makes it a bracket.
205 while (fabs(x1 - x0) >= eps) {
206 double xnew = regulaFalsiStep(x0, x1, f0, f1);
207 double fnew = f(xnew);
208// cout << "xnew = " << xnew << " fnew = " << fnew << endl;
209 ++iterations_used;
210 if (iterations_used > max_iter) {
211 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
212 }
213 if (fabs(fnew) < epsF) {
214 return xnew;
215 }
216 // Now we must check which point we must replace.
217 if ((fnew > 0.0) == (f0 > 0.0)) {
218 // We must replace x0.
219 x0 = x1;
220 f0 = f1;
221 } else {
222 // We must replace x1, this is the case where
223 // the modification to regula falsi kicks in,
224 // by modifying f0.
225 // 1. The classic Illinois method
226// const double gamma = 0.5;
227 // @afr: The next two methods do not work??!!?
228 // 2. The method called 'Method 3' in the paper.
229// const double phi0 = f1/f0;
230// const double phi1 = fnew/f1;
231// const double gamma = 1.0 - phi1/(1.0 - phi0);
232 // 3. The method called 'Method 4' in the paper.
233// const double phi0 = f1/f0;
234// const double phi1 = fnew/f1;
235// const double gamma = 1.0 - phi0 - phi1;
236// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
237// " gamma = " << gamma << endl;
238 // 4. The 'Pegasus' method
239 const double gamma = f1/(f1 + fnew);
240 f0 *= gamma;
241 }
242 x1 = xnew;
243 f1 = fnew;
244 }
245 return 0.5*(x0 + x1);
246 }
247
248
249 private:
250 inline static double regulaFalsiStep(const double a,
251 const double b,
252 const double fa,
253 const double fb)
254 {
255 assert(fa*fb < 0.0);
256 return (b*fa - a*fb)/(fa - fb);
257 }
258
259
260 };
261
262
263
264 template <class ErrorPolicy = ThrowOnError>
266 {
267 public:
268 inline static std::string name()
269 {
270 return "RegulaFalsiBisection";
271 }
272
279 template <class Functor>
280 inline static double solve(const Functor& f,
281 const double a,
282 const double b,
283 const int max_iter,
284 const double tolerance,
285 int& iterations_used)
286 {
287 using namespace std;
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);
291
292 double x0 = a;
293 double x1 = b;
294 double f0 = f(x0);
295 const double epsF = tolerance + macheps * max(fabs(f0), 1.0);
296 if (fabs(f0) < epsF) {
297 return x0;
298 }
299 double f1 = f(x1);
300 if (fabs(f1) < epsF) {
301 return x1;
302 }
303 if (f0 * f1 > 0.0) {
304 return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
305 }
306 iterations_used = 0;
307 // In every iteraton, x1 is the last point computed,
308 // and x0 is the last point computed that makes it a bracket.
309 double width = fabs(x1 - x0);
310 double contraction = 1.0;
311 while (width >= eps) {
312 // If we are contracting sufficiently to at least halve
313 // the interval width in two iterations we use regula
314 // falsi. Otherwise, we take a bisection step to avoid
315 // slow convergence.
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);
320 ++iterations_used;
321 if (iterations_used > max_iter) {
322 return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
323 }
324 if (fabs(fnew) < epsF) {
325 return xnew;
326 }
327 // Now we must check which point we must replace.
328 if ((fnew > 0.0) == (f0 > 0.0)) {
329 // We must replace x0.
330 x0 = x1;
331 f0 = f1;
332 } else {
333 // We must replace x1, this is the case where
334 // the modification to regula falsi kicks in,
335 // by modifying f0.
336 // 1. The classic Illinois method
337 // const double gamma = 0.5;
338 // @afr: The next two methods do not work??!!?
339 // 2. The method called 'Method 3' in the paper.
340 // const double phi0 = f1/f0;
341 // const double phi1 = fnew/f1;
342 // const double gamma = 1.0 - phi1/(1.0 - phi0);
343 // 3. The method called 'Method 4' in the paper.
344 // const double phi0 = f1/f0;
345 // const double phi1 = fnew/f1;
346 // const double gamma = 1.0 - phi0 - phi1;
347 // cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
348 // " gamma = " << gamma << endl;
349 // 4. The 'Pegasus' method
350 const double gamma = f1 / (f1 + fnew);
351 // 5. Straight unmodified Regula Falsi
352 // const double gamma = 1.0;
353 f0 *= gamma;
354 }
355 x1 = xnew;
356 f1 = fnew;
357 const double widthnew = fabs(x1 - x0);
358 contraction = widthnew/width;
359 width = widthnew;
360 }
361 return 0.5 * (x0 + x1);
362 }
363
364
365 private:
366 inline static double regulaFalsiStep(const double a, const double b, const double fa, const double fb)
367 {
368 assert(fa * fb < 0.0);
369 return (b * fa - a * fb) / (fa - fb);
370 }
371 inline static double bisectionStep(const double a, const double b, const double fa, const double fb)
372 {
373 static_cast<void>(fa);
374 static_cast<void>(fb);
375 assert(fa * fb < 0.0);
376 return 0.5*(a + b);
377 }
378 };
379
380
381
382
383
386 template <class Functor>
387 inline void bracketZero(const Functor& f,
388 const double x0,
389 const double dx,
390 double& a,
391 double& b)
392 {
393 const int max_iters = 100;
394 double f0 = f(x0);
395 double cur_dx = dx;
396 int i = 0;
397 for (; i < max_iters; ++i) {
398 double x = x0 + cur_dx;
399 double f_new = f(x);
400 if (f0*f_new <= 0.0) {
401 break;
402 }
403 cur_dx = -2.0*cur_dx;
404 }
405 if (i == max_iters) {
406 OPM_THROW(std::runtime_error,
407 "Could not bracket zero in " +
408 std::to_string(max_iters) + " iterations.");
409 }
410 if (cur_dx < 0.0) {
411 a = x0 + cur_dx;
412 b = i < 2 ? x0 : x0 + 0.25*cur_dx;
413 } else {
414 a = i < 2 ? x0 : x0 + 0.25*cur_dx;
415 b = x0 + cur_dx;
416 }
417 }
418
419} // namespace Opm
420
421#endif // OPM_ROOTFINDERS_HEADER
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