Program Listing for File integration.h

↰ Return to documentation for file (source/utilities/integration.h)

/*
 * integration.h Convenient wrappers around tbs1980's Numerical Integration module for Eigen
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_INTEGRATION_H
#define ES_FLOW_INTEGRATION_H

#include <Eigen/Dense>
#include <Eigen/Core>
#include "NumericalIntegration.h"


namespace utilities {


template<typename T_functor>
Eigen::ArrayXd cumulative_integrate(const Eigen::ArrayXd &x, const T_functor &functor, const int max_subintervals=200)
{
    typedef double ScalarType;
    Eigen::ArrayXd integral(x.rows());
    integral.setZero();

    // Define the integrator and quadrature rule
    Eigen::Integrator<ScalarType> eigIntgtor(max_subintervals);
    Eigen::Integrator<ScalarType>::QuadratureRule quadratureRule = Eigen::Integrator<ScalarType>::GaussKronrod61;

    // Define the desired absolute and relative errors
    auto desAbsErr = ScalarType(0.0);
    ScalarType desRelErr = Eigen::NumTraits<ScalarType>::epsilon() * 50.;

    // Integrate to each value in eta
    for (Eigen::Index i = 0; i < x.rows()-1; i++) {
        integral[i+1] = eigIntgtor.quadratureAdaptive(functor, ScalarType(x[i]), ScalarType(x[i+1]), desAbsErr, desRelErr, quadratureRule) + integral[i];
    }
    return integral;
}


}  /* namespace utilities */

#endif //ES_FLOW_INTEGRATION_H