OGS 6.1.0-1699-ge946d4c5f
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 }
42 
43 
46 template<typename SubType, typename Function>
48 {
49 public:
52  RegulaFalsi(Function const& f, double a, double b)
53  : _f(f), _a(a), _b(b), _fa(f(a)), _fb(f(b))
54  {
55  static_assert(
56  std::is_same<double, decltype(f(0.0))>::value,
57  "Using this class for functions that do not return double"
58  " involves a lot of casts. Hence it is disabled.");
59 
60  if (detail::almost_zero(_fa)) {
61  _b = _a;
62  } else if (detail::almost_zero(_fb)) {
63  _a = _b;
64  } else if (detail::same_sign(_fa, _fb)) {
65  OGS_FATAL("Regula falsi cannot be done, because the function values"
66  " at the interval ends have the same sign.");
67  }
68  }
69 
71  void step(const unsigned num_steps)
72  {
73  for (unsigned i=0; i<num_steps; ++i)
74  {
75  if (_a == _b) return;
76 
77  const double s = (_fb - _fa)/(_b - _a);
78  const double c = _a - _fa/s;
79  const double fc = _f(c);
80 
81  if (detail::almost_zero(fc)) {
82  _a = _b = c;
83  return;
84  }
85  if (!detail::same_sign(fc, _fb))
86  {
87  _a = _b; _fa = _fb;
88  _b = c; _fb = fc;
89  } else {
90  const double m = SubType::get_m(_fa, _fb, fc);
91  _fa *= m;
92  _b = c;
93  _fb = fc;
94  }
95  }
96  }
97 
99  double getResult() const
100  {
101  if (_a == _b) return _a;
102 
103  const double s = (_fb - _fa)/(_b - _a);
104  const double c = _a - _fa/s;
105 
106  return c;
107  }
108 
110  double getRange() const { return std::fabs(_a - _b); }
111 
112 private:
113  Function const& _f;
114  double _a, _b, _fa, _fb;
115 };
116 
117 
123 template<typename SubType, typename Function>
125 makeRegulaFalsi(Function const& f, double const a, double const b)
126 {
127  return RegulaFalsi<SubType, Function>(f, a, b);
128 }
129 
130 
133 {
134  static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
135  { return 1.0; }
136 };
137 
139 struct Illinois
140 {
141  static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
142  { return 0.5; }
143 };
144 
146 struct Pegasus
147 {
148  static double get_m(const double /*fa*/, const double fb, const double fc)
149  { return fb / (fb+fc); }
150 };
151 
154 {
155  static double get_m(const double /*fa*/, const double fb, const double fc)
156  {
157  const double f = 1.0 - fc / fb;
158  return (f >= 0.0) ? f : 0.5;
159  }
160 };
161 
162 }
163 
164 } // namespace MathLib
double getRange() const
Returns the size of the current search interval.
Definition: Root1D.h:110
bool almost_zero(double a)
Definition: Root1D.h:36
Used by RegulaFalsi in the original regula falsi algorithm.
Definition: Root1D.h:132
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:134
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:139
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:146
RegulaFalsi(Function const &f, double a, double b)
Definition: Root1D.h:52
double getResult() const
Returns the current estimate of the root.
Definition: Root1D.h:99
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:155
Used by RegulaFalsi in a modified version of the regula falsi algorithm.
Definition: Root1D.h:153
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void step(const unsigned num_steps)
Do num_steps iteration of regula falsi.
Definition: Root1D.h:71
static double get_m(const double, const double, const double)
Definition: Root1D.h:141
RegulaFalsi< SubType, Function > makeRegulaFalsi(Function const &f, double const a, double const b)
Definition: Root1D.h:125
static double get_m(const double, const double fb, const double fc)
Definition: Root1D.h:148