OGS 6.2.1-97-g73d1aeda3
Root1D.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <cassert>
13 #include <cmath>
14 #include <limits>
15 #include <type_traits>
16 #include <logog/include/logog.hpp>
17 #include "BaseLib/Error.h"
18 
19 namespace MathLib
20 {
21 namespace Nonlinear
22 {
23 
24 namespace detail
25 {
26 
28 inline bool same_sign(double a, double b)
29 {
30  // note: signbit() casts integers to double and distinguishes between +0.0 and
31  // -0.0. However, the latter is not a problem, because if f(x0) == 0, we've
32  // already found the root x0 that we are searching for.
33  return std::signbit(a) == std::signbit(b);
34 }
35 
36 inline bool almost_zero(double a)
37 {
38  return std::abs(a) <= std::numeric_limits<double>::epsilon();
39 }
40 
41 } // namespace detail
42 
45 template<typename SubType, typename Function>
47 {
48 public:
51  RegulaFalsi(Function const& f, double a, double b)
52  : _f(f), _a(a), _b(b), _fa(f(a)), _fb(f(b))
53  {
54  static_assert(
55  std::is_same<double, decltype(f(0.0))>::value,
56  "Using this class for functions that do not return double"
57  " involves a lot of casts. Hence it is disabled.");
58 
59  if (detail::almost_zero(_fa)) {
60  _b = _a;
61  } else if (detail::almost_zero(_fb)) {
62  _a = _b;
63  } else if (detail::same_sign(_fa, _fb)) {
64  OGS_FATAL("Regula falsi cannot be done, because the function values"
65  " at the interval ends have the same sign.");
66  }
67  }
68 
70  void step(const unsigned num_steps)
71  {
72  for (unsigned i=0; i<num_steps; ++i)
73  {
74  if (_a == _b)
75  {
76  return;
77  }
78 
79  const double s = (_fb - _fa)/(_b - _a);
80  const double c = _a - _fa/s;
81  const double fc = _f(c);
82 
83  if (detail::almost_zero(fc)) {
84  _a = _b = c;
85  return;
86  }
87  if (!detail::same_sign(fc, _fb))
88  {
89  _a = _b; _fa = _fb;
90  _b = c; _fb = fc;
91  } else {
92  const double m = SubType::get_m(_fa, _fb, fc);
93  _fa *= m;
94  _b = c;
95  _fb = fc;
96  }
97  }
98  }
99 
101  double getResult() const
102  {
103  if (_a == _b)
104  {
105  return _a;
106  }
107 
108  const double s = (_fb - _fa)/(_b - _a);
109  const double c = _a - _fa/s;
110 
111  return c;
112  }
113 
115  double getRange() const { return std::fabs(_a - _b); }
116 
117 private:
118  Function const& _f;
119  double _a, _b, _fa, _fb;
120 };
121 
122 
128 template<typename SubType, typename Function>
130 makeRegulaFalsi(Function const& f, double const a, double const b)
131 {
132  return RegulaFalsi<SubType, Function>(f, a, b);
133 }
134 
135 
138 {
139  static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
140  { return 1.0; }
141 };
142 
144 struct Illinois
145 {
146  static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
147  { return 0.5; }
148 };
149 
151 struct Pegasus
152 {
153  static double get_m(const double /*fa*/, const double fb, const double fc)
154  { return fb / (fb+fc); }
155 };
156 
159 {
160  static double get_m(const double /*fa*/, const double fb, const double fc)
161  {
162  const double f = 1.0 - fc / fb;
163  return (f >= 0.0) ? f : 0.5;
164  }
165 };
166 
167 } // namespace Nonlinear
168 
169 } // namespace MathLib
double getRange() const
Returns the size of the current search interval.
Definition: Root1D.h:115
bool almost_zero(double a)
Definition: Root1D.h:36
Used by RegulaFalsi in the original regula falsi algorithm.
Definition: Root1D.h:137
std::function< bool(const double t, MappedConstVector< N > const &y, MappedVector< N > &ydot)> Function
static double get_m(const double, const double, const double)
Definition: Root1D.h:139
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:144
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:151
RegulaFalsi(Function const &f, double a, double b)
Definition: Root1D.h:51
double getResult() const
Returns the current estimate of the root.
Definition: Root1D.h:101
bool same_sign(double a, double b)
Tells if a and b have the same sign.
Definition: Root1D.h:28
static double get_m(const double, const double fb, const double fc)
Definition: Root1D.h:160
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:158
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void step(const unsigned num_steps)
Do num_steps iteration of regula falsi.
Definition: Root1D.h:70
static double get_m(const double, const double, const double)
Definition: Root1D.h:146
RegulaFalsi< SubType, Function > makeRegulaFalsi(Function const &f, double const a, double const b)
Definition: Root1D.h:130
static double get_m(const double, const double fb, const double fc)
Definition: Root1D.h:153