NumCpp  2.10.1
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::roots
43 {
44  //================================================================================
45  // Class Description:
48  class Dekker : public Iteration
49  {
50  public:
51  //============================================================================
52  // Method Description:
58  Dekker(const double epsilon, std::function<double(double)> f) noexcept :
59  Iteration(epsilon),
60  f_(std::move(f))
61  {
62  }
63 
64  //============================================================================
65  // Method Description:
72  Dekker(const double epsilon, const uint32 maxNumIterations, std::function<double(double)> f) noexcept :
73  Iteration(epsilon, maxNumIterations),
74  f_(std::move(f))
75  {
76  }
77 
78  //============================================================================
79  // Method Description:
82  ~Dekker() override = default;
83 
84  //============================================================================
85  // Method Description:
92  double solve(double a, double b)
93  {
95 
96  double fa = f_(a);
97  double fb = f_(b);
98 
99  checkAndFixAlgorithmCriteria(a, b, fa, fb);
100 
101  double lastB = a;
102  double lastFb = fa;
103 
104  while (std::fabs(fb) > epsilon_ && std::fabs(b - a) > epsilon_)
105  {
106  const double s = calculateSecant(b, fb, lastB, lastFb);
107  const double m = calculateBisection(a, b);
108 
109  lastB = b;
110 
111  b = useSecantMethod(b, s, m) ? s : m;
112 
113  lastFb = fb;
114  fb = f_(b);
115 
116  if (fa * fb > 0 && fb * lastFb < 0)
117  {
118  a = lastB;
119  }
120 
121  fa = f_(a);
122  checkAndFixAlgorithmCriteria(a, b, fa, fb);
123 
125  }
126 
127  return b;
128  }
129 
130  private:
131  //============================================================================
132  const std::function<double(double)> f_;
133 
134  //============================================================================
135  // Method Description:
143  static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) noexcept
144  {
145  // Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
146  if (std::fabs(fa) < std::fabs(fb))
147  {
148  std::swap(a, b);
149  std::swap(fa, fb);
150  }
151  }
152 
153  //============================================================================
154  // Method Description:
163  static double calculateSecant(double b, double fb, double lastB, double lastFb) noexcept
164  {
165  // No need to check division by 0, in this case the method returns NAN which is taken care by
166  // useSecantMethod method
167  return b - fb * (b - lastB) / (fb - lastFb);
168  }
169 
170  //============================================================================
171  // Method Description:
178  static double calculateBisection(double a, double b) noexcept
179  {
180  return 0.5 * (a + b);
181  }
182 
183  //============================================================================
184  // Method Description:
192  static bool useSecantMethod(double b, double s, double m) noexcept
193  {
194  // Value s calculated by secant method has to be between m and b
195  return (b > m && s > m && s < b) || (b < m && s > b && s < m);
196  }
197  };
198 } // namespace nc::roots
Definition: Dekker.hpp:49
~Dekker() override=default
double solve(double a, double b)
Definition: Dekker.hpp:92
Dekker(const double epsilon, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:58
Dekker(const double epsilon, const uint32 maxNumIterations, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:72
ABC for iteration classes to derive from.
Definition: Iteration.hpp:46
Iteration(double epsilon) noexcept
Definition: Iteration.hpp:54
const double epsilon_
Definition: Iteration.hpp:116
void resetNumberOfIterations() noexcept
Definition: Iteration.hpp:94
void incrementNumberOfIterations()
Definition: Iteration.hpp:105
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:56
Definition: Bisection.hpp:43
void swap(NdArray< dtype > &inArray1, NdArray< dtype > &inArray2) noexcept
Definition: swap.hpp:42
std::uint32_t uint32
Definition: Types.hpp:40