Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
std::vector<std::vector<double>> rombergIntegral(double a, double b, size_t n, const std::function<double (double)> &f) {
std::vector<std::vector<double>> romberg_integral(n, std::vector<double>(n));
//R(0,0) Start with trapezoidal integration with N = 1
romberg_integral.front().front() = trapezoidalIntegral(a, b, 1, f);
double h = b-a;
for(size_t step = 1; step < n; step++) {
h *= 0.5;
//R(step, 0) Improve trapezoidal integration with decreasing h
double trapezoidal_integration = 0;
size_t stepEnd = pow(2, step - 1);
for(size_t tzStep = 1; tzStep <= stepEnd; tzStep++) {
const double deltaX = (2*tzStep - 1)*h;
trapezoidal_integration += f(a + deltaX);
}
romberg_integral[step].front() = 0.5*romberg_integral[step - 1].front() + trapezoidal_integration*h;
//R(m,n) Romberg integration with R(m,1) -> Simpson rule, R(m,2) -> Boole's rule
for(size_t rbStep = 1; rbStep <= step; rbStep++) {
const double k = pow(4, rbStep);
romberg_integral[step][rbStep] = (k*romberg_integral[step][rbStep-1] - romberg_integral[step-1][rbStep-1])/(k-1);
}
}
return romberg_integral;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.