adaptive_simpson.hpp¶
Numerically integrate a smooth real function on [left, right] with adaptive Simpson subdivision and an absolute error target.
按误差估计自适应细分区间，数值计算平滑函数的定积分。
Implementation¶
#ifndef NOYA_ADAPTIVE_SIMPSON_HPP
#define NOYA_ADAPTIVE_SIMPSON_HPP 1
/// @complexity Time: O(2^d) function evaluations in the worst case, where d is max_depth.
/// Space: O(d) recursion stack.
#include <cassert>
#include <cmath>
#include <functional>
namespace noya {
/// @brief Numerically integrate a smooth real function on [left, right] with
/// adaptive Simpson subdivision and an absolute error target.
template <class Function, class Real = long double>
Real adaptive_simpson(Function function, Real left, Real right,
Real epsilon = Real(1e-10), int max_depth = 24) {
assert(epsilon > Real{});
assert(max_depth >= 0);
if (left == right) {
return Real{};
}
if (right < left) {
return -adaptive_simpson<Function, Real>(function, right, left, epsilon,
max_depth);
}
auto simpson = [](Real first, Real second, Real f_first, Real f_middle,
Real f_second) {
return (second - first) * (f_first + Real(4) * f_middle + f_second) /
Real(6);
};
Real middle = (left + right) / Real(2);
Real f_left = std::invoke(function, left);
Real f_middle = std::invoke(function, middle);
Real f_right = std::invoke(function, right);
Real whole = simpson(left, right, f_left, f_middle, f_right);
auto solve = [&](auto &self, Real first, Real second, Real f_first,
Real f_center, Real f_second, Real estimate, Real tolerance,
int depth) -> Real {
Real center = (first + second) / Real(2);
Real left_center = (first + center) / Real(2);
Real right_center = (center + second) / Real(2);
Real f_left_center = std::invoke(function, left_center);
Real f_right_center = std::invoke(function, right_center);
Real left_estimate =
simpson(first, center, f_first, f_left_center, f_center);
Real right_estimate =
simpson(center, second, f_center, f_right_center, f_second);
Real difference = left_estimate + right_estimate - estimate;
if (depth == 0 || std::abs(difference) <= Real(15) * tolerance) {
return left_estimate + right_estimate + difference / Real(15);
}
return self(self, first, center, f_first, f_left_center, f_center,
left_estimate, tolerance / Real(2), depth - 1) +
self(self, center, second, f_center, f_right_center, f_second,
right_estimate, tolerance / Real(2), depth - 1);
};
return solve(solve, left, right, f_left, f_middle, f_right, whole, epsilon,
max_depth);
}
} // namespace noya
#endif // NOYA_ADAPTIVE_SIMPSON_HPP
#include <cassert>
#include <cmath>
#include <functional>
/// @complexity Time: O(2^d) function evaluations in the worst case, where d is max_depth.
/// Space: O(d) recursion stack.
namespace noya {
/// @brief Numerically integrate a smooth real function on [left, right] with
/// adaptive Simpson subdivision and an absolute error target.
template <class Function, class Real = long double>
Real adaptive_simpson(Function function, Real left, Real right,
Real epsilon = Real(1e-10), int max_depth = 24) {
assert(epsilon > Real{});
assert(max_depth >= 0);
if (left == right) {
return Real{};
}
if (right < left) {
return -adaptive_simpson<Function, Real>(function, right, left, epsilon,
max_depth);
}
auto simpson = [](Real first, Real second, Real f_first, Real f_middle,
Real f_second) {
return (second - first) * (f_first + Real(4) * f_middle + f_second) /
Real(6);
};
Real middle = (left + right) / Real(2);
Real f_left = std::invoke(function, left);
Real f_middle = std::invoke(function, middle);
Real f_right = std::invoke(function, right);
Real whole = simpson(left, right, f_left, f_middle, f_right);
auto solve = [&](auto &self, Real first, Real second, Real f_first,
Real f_center, Real f_second, Real estimate, Real tolerance,
int depth) -> Real {
Real center = (first + second) / Real(2);
Real left_center = (first + center) / Real(2);
Real right_center = (center + second) / Real(2);
Real f_left_center = std::invoke(function, left_center);
Real f_right_center = std::invoke(function, right_center);
Real left_estimate =
simpson(first, center, f_first, f_left_center, f_center);
Real right_estimate =
simpson(center, second, f_center, f_right_center, f_second);
Real difference = left_estimate + right_estimate - estimate;
if (depth == 0 || std::abs(difference) <= Real(15) * tolerance) {
return left_estimate + right_estimate + difference / Real(15);
}
return self(self, first, center, f_first, f_left_center, f_center,
left_estimate, tolerance / Real(2), depth - 1) +
self(self, center, second, f_center, f_right_center, f_second,
right_estimate, tolerance / Real(2), depth - 1);
};
return solve(solve, left, right, f_left, f_middle, f_right, whole, epsilon,
max_depth);
}
} // namespace noya