NumCpp  2.1.0
A C++ implementation of the Python Numpy library
romberg.hpp
Go to the documentation of this file.
1 #pragma once
34 
35 #include "NumCpp/Core/Types.hpp"
37 #include "NumCpp/Utils/power.hpp"
38 
39 #include <functional>
40 #include <vector>
41 
42 namespace nc
43 {
44  namespace integrate
45  {
46  //============================================================================
47  // Method Description:
57  inline double romberg(const double low, const double high, const uint8 n,
58  const std::function<double(double)>& f)
59  {
60  NdArray<double> rombergIntegral(n);
61 
62  //R(0,0) Start with trapezoidal integration with N = 1
63  rombergIntegral(0, 0) = trapazoidal(low, high, 1, f);
64 
65  double h = high - low;
66  for (uint8 step = 1; step < n; step++)
67  {
68  h *= 0.5;
69 
70  //R(step, 0) Improve trapezoidal integration with decreasing h
71  double trapezoidal_integration = 0.0;
72  const uint32 stepEnd = utils::power(2, step - 1);
73  for (uint32 tzStep = 1; tzStep <= stepEnd; ++tzStep)
74  {
75  const double deltaX = (2.0 * static_cast<double>(tzStep - 1)) * h;
76  trapezoidal_integration += f(low + deltaX);
77  }
78 
79  rombergIntegral(step, 0) = 0.5 * rombergIntegral(step - 1, 0);
80  rombergIntegral(step, 0) += trapezoidal_integration * h;
81 
82  //R(m,n) Romberg integration with R(m,1) -> Simpson rule, R(m,2) -> Boole's rule
83  for (uint8 rbStep = 1; rbStep <= step; ++rbStep)
84  {
85  const double k = utils::power(4, rbStep);
86  rombergIntegral(step, rbStep) = k * rombergIntegral(step, rbStep - 1);
87  rombergIntegral(step, rbStep) -= rombergIntegral(step - 1, rbStep - 1);
88  rombergIntegral(step, rbStep) /= (k - 1.0);
89  }
90  }
91 
92  return rombergIntegral.back();
93  }
94  } // namespace integrate
95 } // namespace nc
nc::uint8
std::uint8_t uint8
Definition: Types.hpp:43
nc::NdArray< double >
nc::NdArray::back
value_type back() const noexcept
Definition: NdArrayCore.hpp:2221
nc::uint32
std::uint32_t uint32
Definition: Types.hpp:41
trapazoidal.hpp
power.hpp
nc
Definition: Coordinate.hpp:45
nc::utils::power
dtype power(dtype inValue, uint8 inPower) noexcept
Definition: Utils/power.hpp:49
nc::integrate::trapazoidal
double trapazoidal(const double low, const double high, const uint32 n, const std::function< double(double)> &f) noexcept
Definition: trapazoidal.hpp:54
Types.hpp
nc::random::f
dtype f(dtype inDofN, dtype inDofD)
Definition: f.hpp:58
nc::integrate::romberg
double romberg(const double low, const double high, const uint8 n, const std::function< double(double)> &f)
Definition: romberg.hpp:57