NumCpp  2.1.0
A C++ implementation of the Python Numpy library
Dekker.hpp
Go to the documentation of this file.
1 #pragma once
34 
35 #include "NumCpp/Core/Types.hpp"
37 
38 #include <cmath>
39 #include <functional>
40 #include <utility>
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,
61  std::function<double(double)> f) noexcept :
62  Iteration(epsilon),
63  f_(std::move(f))
64  {}
65 
66  //============================================================================
67  // Method Description:
74  Dekker(const double epsilon,
75  const uint32 maxNumIterations,
76  std::function<double(double)> f) noexcept :
77  Iteration(epsilon, maxNumIterations),
78  f_(std::move(f))
79  {}
80 
81  //============================================================================
82  // Method Description:
85  ~Dekker() override = default;
86 
87  //============================================================================
88  // Method Description:
95  double solve(double a, double b)
96  {
98 
99  double fa = f_(a);
100  double fb = f_(b);
101 
102  checkAndFixAlgorithmCriteria(a, b, fa, fb);
103 
104  double lastB = a;
105  double lastFb = fa;
106 
107  while (std::fabs(fb) > epsilon_ && std::fabs(b - a) > epsilon_)
108  {
109  const double s = calculateSecant(b, fb, lastB, lastFb);
110  const double m = calculateBisection(a, b);
111 
112  lastB = b;
113 
114  b = useSecantMethod(b, s, m) ? s : m;
115 
116  lastFb = fb;
117  fb = f_(b);
118 
119  if (fa * fb > 0 && fb * lastFb < 0)
120  {
121  a = lastB;
122  }
123 
124  fa = f_(a);
125  checkAndFixAlgorithmCriteria(a, b, fa, fb);
126 
128  }
129 
130  return b;
131  }
132 
133  private:
134  //============================================================================
135  const std::function<double(double)> f_;
136 
137  //============================================================================
138  // Method Description:
146  static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) noexcept
147  {
148  //Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
149  if (std::fabs(fa) < std::fabs(fb))
150  {
151  std::swap(a, b);
152  std::swap(fa, fb);
153  }
154  }
155 
156  //============================================================================
157  // Method Description:
166  static double calculateSecant(double b, double fb, double lastB, double lastFb) noexcept
167  {
168  //No need to check division by 0, in this case the method returns NAN which is taken care by 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) ||
198  (b < m && s > b && s < m);
199  }
200  };
201  } // namespace roots
202 } // namespace nc
nc::roots::Dekker::Dekker
Dekker(const double epsilon, const uint32 maxNumIterations, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:74
nc::roots::Iteration::Iteration
Iteration(double epsilon) noexcept
Definition: Iteration.hpp:56
nc::roots::Dekker::solve
double solve(double a, double b)
Definition: Dekker.hpp:95
nc::roots::Dekker
Definition: Dekker.hpp:50
nc::roots::Iteration
ABC for iteration classes to derive from.
Definition: Iteration.hpp:47
nc::roots::Iteration::incrementNumberOfIterations
void incrementNumberOfIterations()
Definition: Iteration.hpp:105
nc::uint32
std::uint32_t uint32
Definition: Types.hpp:41
nc::roots::Dekker::Dekker
Dekker(const double epsilon, std::function< double(double)> f) noexcept
Definition: Dekker.hpp:60
nc::roots::Iteration::epsilon_
const double epsilon_
Definition: Iteration.hpp:115
Iteration.hpp
nc::roots::Iteration::resetNumberOfIterations
void resetNumberOfIterations() noexcept
Definition: Iteration.hpp:94
nc
Definition: Coordinate.hpp:45
nc::swap
void swap(NdArray< dtype > &inArray1, NdArray< dtype > &inArray2) noexcept
Definition: swap.hpp:43
nc::roots::Dekker::~Dekker
~Dekker() override=default
Types.hpp
nc::random::f
dtype f(dtype inDofN, dtype inDofD)
Definition: f.hpp:58