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