NumCpp  2.8.0
A Templatized Header Only C++ Implementation of the Python NumPy Library
Dekker.hpp
Go to the documentation of this file.
1 #pragma once
34 
35 #include <cmath>
36 #include <functional>
37 #include <utility>
38 
39 #include "NumCpp/Core/Types.hpp"
41 
42 namespace nc
43 {
44  namespace roots
45  {
46  //================================================================================
47  // Class Description:
50  class Dekker : public Iteration
51  {
52  public:
53  //============================================================================
54  // Method Description:
60  Dekker(const double epsilon, std::function<double(double)> f) noexcept :
61  Iteration(epsilon),
62  f_(std::move(f))
63  {
64  }
65 
66  //============================================================================
67  // Method Description:
74  Dekker(const double epsilon, const uint32 maxNumIterations, std::function<double(double)> f) noexcept :
75  Iteration(epsilon, maxNumIterations),
76  f_(std::move(f))
77  {
78  }
79 
80  //============================================================================
81  // Method Description:
84  ~Dekker() override = default;
85 
86  //============================================================================
87  // Method Description:
94  double solve(double a, double b)
95  {
97 
98  double fa = f_(a);
99  double fb = f_(b);
100 
101  checkAndFixAlgorithmCriteria(a, b, fa, fb);
102 
103  double lastB = a;
104  double lastFb = fa;
105 
106  while (std::fabs(fb) > epsilon_ && std::fabs(b - a) > epsilon_)
107  {
108  const double s = calculateSecant(b, fb, lastB, lastFb);
109  const double m = calculateBisection(a, b);
110 
111  lastB = b;
112 
113  b = useSecantMethod(b, s, m) ? s : m;
114 
115  lastFb = fb;
116  fb = f_(b);
117 
118  if (fa * fb > 0 && fb * lastFb < 0)
119  {
120  a = lastB;
121  }
122 
123  fa = f_(a);
124  checkAndFixAlgorithmCriteria(a, b, fa, fb);
125 
127  }
128 
129  return b;
130  }
131 
132  private:
133  //============================================================================
134  const std::function<double(double)> f_;
135 
136  //============================================================================
137  // Method Description:
145  static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) noexcept
146  {
147  // Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
148  if (std::fabs(fa) < std::fabs(fb))
149  {
150  std::swap(a, b);
151  std::swap(fa, fb);
152  }
153  }
154 
155  //============================================================================
156  // Method Description:
165  static double calculateSecant(double b, double fb, double lastB, double lastFb) noexcept
166  {
167  // No need to check division by 0, in this case the method returns NAN which is taken care by
168  // useSecantMethod method
169  return b - fb * (b - lastB) / (fb - lastFb);
170  }
171 
172  //============================================================================
173  // Method Description:
180  static double calculateBisection(double a, double b) noexcept
181  {
182  return 0.5 * (a + b);
183  }
184 
185  //============================================================================
186  // Method Description:
194  static bool useSecantMethod(double b, double s, double m) noexcept
195  {
196  // Value s calculated by secant method has to be between m and b
197  return (b > m && s > m && s < b) || (b < m && s > b && s < m);
198  }
199  };
200  } // namespace roots
201 } // namespace nc
Definition: Dekker.hpp:51
~Dekker() override=default
double solve(double a, double b)
Definition: Dekker.hpp:94
Dekker(const double epsilon, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:60
Dekker(const double epsilon, const uint32 maxNumIterations, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:74
ABC for iteration classes to derive from.
Definition: Iteration.hpp:48
Iteration(double epsilon) noexcept
Definition: Iteration.hpp:56
const double epsilon_
Definition: Iteration.hpp:118
void resetNumberOfIterations() noexcept
Definition: Iteration.hpp:96
void incrementNumberOfIterations()
Definition: Iteration.hpp:107
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:58
Definition: Coordinate.hpp:45
void swap(NdArray< dtype > &inArray1, NdArray< dtype > &inArray2) noexcept
Definition: swap.hpp:42
std::uint32_t uint32
Definition: Types.hpp:40