37 double integrate(
const std::function<
double(
double)>& _f,
const double _a,
const double _b,
double _eps,
const uint _minIter,
const uint _maxIter,
const uint _branchFactor,
const uint _maxRecursion,
const bool _relativeError) {
40 double min = std::min(_a,_b);
41 double max = std::max(_a,_b);
42 if (_relativeError) { _eps = std::max(_eps, std::numeric_limits<double>::epsilon()); }
43 std::vector<double> iterants;
44 iterants.push_back((max - min)*(_f(min) + _f(max))/2);
48 double maxVal = std::abs(iterants[0]);
49 for (
uint iter=0; iter<_maxIter; ++iter) {
51 for (
double x = min+h/2; x<max; x += h) {
54 maxVal = std::max(maxVal, std::abs(f));
59 sum += iterants.back()/2;
60 iterants.push_back(sum);
61 double oldIt0 = iterants[0];
62 for (
size_t k=0; k<iterants.size()-1; ++k) {
63 size_t i=iterants.size()-1-k;
64 iterants[i-1] = iterants[i] + (iterants[i]-iterants[i-1])/(
misc::pow(2.0, 2*(k+1))-1);
67 error = std::abs((iterants[0]-oldIt0)/oldIt0);
68 if (std::isnan(error)) { error = std::abs(iterants[0]-oldIt0); }
70 error = std::abs(iterants[0]-oldIt0);
72 if (iter >= _minIter && error <= _eps) {
73 return (_a>_b?-1:1)*iterants[0];
76 if (_maxRecursion == 0) {
77 return (_a>_b?-1:1)*iterants[0];
81 LOG(debug, (_relativeError?
"relative":
"absolute") <<
" error " << error <<
" is larger than eps = " << _eps);
83 h = (max-min) / _branchFactor;
84 double newEps = (_relativeError
85 ? std::max(std::abs(iterants[0]), maxVal)*_eps
88 std::sqrt(_branchFactor)*std::numeric_limits<double>::epsilon()*std::max(std::abs(iterants[0]), maxVal)
91 for (
uint i=0; i<_branchFactor; ++i) {
92 sum +=
integrate(_f, min+i*h, min+(i+1)*h, newEps, _minIter, _maxIter, _branchFactor, _maxRecursion-1,
false);
94 return (_a>_b?-1:1)*sum;
98 double integrate_segmented(
const std::function<
double(
double)> &_f,
const double _a,
const double _b,
const double _segmentation,
const double _eps,
const uint _minIter,
const uint _maxIter,
const uint _branchFactor,
const uint _maxRecursion) {
99 const double min = std::min(_a, _b);
100 const double max = std::max(_a, _b);
102 for (
double x = min; x < max; x += _segmentation) {
103 res +=
integrate(_f, x, std::min(x+_segmentation, max), _eps, _minIter, _maxIter, _branchFactor, _maxRecursion);
105 return (_a>_b?-1:1)*res;
109 double find_root_bisection(
const std::function<
double(
double)> &_f,
double _min,
double _max,
const double _epsilon) {
110 const double nm = std::max(_min, _max);
111 _min = std::min(_min, _max);
114 double fmin = _f(_min);
115 double fmax = _f(_max);
117 REQUIRE(fmin * fmax <= 0,
"bisection requires inputs to both sides of the root (and no NaNs)");
119 if (std::fpclassify(fmin) == FP_ZERO) {
123 if (std::fpclassify(fmax) == FP_ZERO) {
127 while (_max-_min > _epsilon) {
128 double mid = (_max+_min)/2.0;
129 double fmid = _f(mid);
130 if (std::fpclassify(fmid) == FP_ZERO) {
133 REQUIRE(std::isfinite(fmid),
"invalid function value f("<<mid<<
") = " << fmid <<
" reached in bisection");
134 if (fmin * fmid < 0) {
143 return (_max + _min) / 2.0;
148 Polynomial::Polynomial() {}
150 Polynomial::Polynomial(std::vector<double> _coeff) : coefficients(std::move(_coeff)) {}
152 size_t Polynomial::terms()
const {
153 return coefficients.size();
156 Polynomial& Polynomial::operator+=(
const Polynomial &_rhs) {
157 coefficients.resize(std::max(terms(), _rhs.terms()));
159 for(
size_t i = 0; i < _rhs.terms(); ++i) {
160 coefficients[i] += _rhs.coefficients[i];
166 Polynomial& Polynomial::operator-=(
const Polynomial &_rhs) {
167 coefficients.resize(std::max(terms(), _rhs.terms()));
169 for(
size_t i = 0; i < _rhs.terms(); ++i) {
170 coefficients[i] -= _rhs.coefficients[i];
177 Polynomial& Polynomial::operator/=(
double _rhs) {
178 for (
size_t i=0; i<terms(); ++i) {
179 coefficients[i] /= _rhs;
184 Polynomial& Polynomial::operator*=(
double _rhs) {
185 for (
size_t i=0; i<terms(); ++i) {
186 coefficients[i] *= _rhs;
193 result.coefficients.resize(_rhs.terms()+terms()-1);
194 for (
size_t i=0; i<terms(); ++i) {
195 for (
size_t j=0; j<_rhs.terms(); ++j) {
196 result.coefficients[i+j] += coefficients[i]*_rhs.coefficients[j];
204 result.coefficients.resize(terms());
205 for (
size_t i=0; i<terms(); ++i) {
206 result.coefficients[i] = coefficients[i]*_rhs;
211 double Polynomial::operator()(
double x)
const {
213 for (
size_t i=terms(); i>0; --i) {
215 result += coefficients[i-1];
220 double Polynomial::scalar_product(
const Polynomial &_rhs,
const std::function<
double (
double)> &_weight,
double _minX,
double _maxX)
const {
222 return (*
this)(x) * _rhs(x) * _weight(x);
223 }, _minX, _maxX, 1e-10);
226 double Polynomial::norm(
const std::function<
double (
double)> &_weight,
double _minX,
double _maxX)
const {
227 return std::sqrt(scalar_product(*
this, _weight, _minX, _maxX));
231 Polynomial& Polynomial::orthogonolize(
const std::vector<Polynomial> &_orthoBase,
const std::function<
double (
double)> &_weight,
double _minX,
double _maxX) {
232 for (
size_t i=0; i<_orthoBase.size(); ++i) {
233 (*this)-=_orthoBase[i]*scalar_product(_orthoBase[i], _weight, _minX, _maxX);
234 REQUIRE(std::abs(scalar_product(_orthoBase[i], _weight, _minX, _maxX)) < 1e-12, i <<
" " << std::abs(scalar_product(_orthoBase[i], _weight, _minX, _maxX)));
236 (*this)/=norm(_weight, _minX, _maxX);
241 std::vector<Polynomial> Polynomial::build_orthogonal_base(
size_t _N,
const std::function<
double (
double)> &_weight,
double _minX,
double _maxX) {
242 std::vector<Polynomial> base;
243 while (base.size() < _N) {
245 next.coefficients.resize(base.size()+1);
246 next.coefficients.back() = 1;
247 next.orthogonolize(base, _weight, _minX, _maxX);
248 base.push_back(next);
249 LOG(debug,
"next basis polynomial " << next.coefficients);
257 template<
class ft_type>
259 ft_type denominator = x1 - 2*x2 + x3;
260 if (std::abs(denominator) < epsilon * std::max(x1, std::max(x2, x3))) {
263 return (x1*x3 - x2*x2) / denominator;
266 template<
class ft_type>
268 values.push_back(_val);
269 for (
size_t i=values.size()-1; i>=2; i-=2) {
270 values[i-2] = shanks(values[i-2], values[i-1], values[i]);
274 template<
class ft_type>
276 if (values.empty()) {
279 return values[(values.size()-1) % 2];
283 template<
class ft_type>
285 size_t i = (values.size()-1) % 2;
286 if (i+1 >= values.size()) {
289 return std::abs(values[i] - values[i+1]);
292 template<
class ft_type>
305 template<
class ft_type>
307 return ft_type(n+1)*x2 - ft_type(n)*x1;
310 template<
class ft_type>
312 values.push_back(_val);
313 for (
size_t i=values.size()-1; i>=1; i-=1) {
314 values[i-1] = richard(i-1, values[i-1], values[i]);
318 template<
class ft_type>
320 if (values.empty()) {
323 return values.front();
327 template<
class ft_type>
329 if (values.size() < 2) {
332 return std::abs(values[0] - values[1]);
335 template<
class ft_type>
Header file for some additional math functions.
Header file for CHECK and REQUIRE macros.
double integrate_segmented(const std::function< double(double)> &_f, double _a, double _b, double _segmentation, double _eps=1e-8, uint _minIter=4, uint _maxIter=6, uint _branchFactor=8, uint _maxRecursion=10)
Header file for the standard contaienr support functions.
The main namespace of xerus.
Header file for several elementary numerical methods: polynomials, romberg integration and limit extr...
The xerus exception class.
Header file for xerus::misc::generic_error exception class.
Header file for all logging macros and log-buffer functionality.
double find_root_bisection(const std::function< double(double)> &_f, double _min, double _max, double _epsilon=1e-14)
Header file for comfort functions and macros that should not be exported in the library.
#define XERUS_THROW(...)
Helper macro to throw a generic_error (or derived exception) with some additional information include...
double integrate(const std::function< double(double)> &_f, double _a, double _b, double _eps=std::numeric_limits< double >::epsilon(), uint _minIter=4, uint _maxIter=6, uint _branchFactor=7, uint _maxRecursion=10, bool _relativeError=true)
Performs a Romberg Integration (richardson extrapolation of regular riemannian sum) + adaptive refine...
IndexedTensorMoveable< tensor_type > operator*(const value_t _factor, IndexedTensorMoveable< tensor_type > &&_tensor)
constexpr T pow(const T &_base, const uint32 _exp) noexcept
: Calculates _base^_exp by binary exponentiation