Put in its own method to avoid Eigen dependency

This commit is contained in:
Ian Bell
2025-05-20 07:56:06 -04:00
parent 8e6f2ef03f
commit c419eff1ec
2 changed files with 24 additions and 23 deletions

View File

@@ -710,27 +710,4 @@ double asinh(double value) {
#endif
template<typename Callable>
auto romberg_diff(Callable& func, double x, std::size_t order=2, int h=0.1){
// Initialize the table
auto r = Eigen::ArrayXd(order, order);
// Compute the first column using the central difference formula
for (auto i = 0; i < order; ++i){
r(i, 0) = (func(x + h) - func(x - h)) / (2 * h);
h /= 2.0;
}
// Apply Richardson extrapolation
for (auto i = 1; i < order; ++i){
for (auto j = i; j < order; ++j){
double fouri = pow(4, i);
r(j, i) = (fouri * r(j, i-1) - r(j-1, i-1)) / (fouri - 1);
}
}
return r(order - 1, order - 1);
}
#endif

24
include/finite_diff.h Normal file
View File

@@ -0,0 +1,24 @@
#include <Eigen/Dense>
template<typename Callable>
auto romberg_diff(Callable& func, double x, std::size_t order=2, double h=0.1){
// Initialize the table
auto r = Eigen::ArrayXd(order, order);
// Compute the first column using the central difference formula
for (auto i = 0; i < order; ++i){
r(i, 0) = (func(x + h) - func(x - h)) / (2 * h);
h /= 2.0;
}
// Apply Richardson extrapolation
for (auto i = 1; i < order; ++i){
for (auto j = i; j < order; ++j){
double fouri = pow(4, i);
r(j, i) = (fouri * r(j, i-1) - r(j-1, i-1)) / (fouri - 1);
}
}
return r(order - 1, order - 1);
}