mirror of
https://github.com/gwm17/catima.git
synced 2024-11-22 18:28:51 -05:00
105 lines
2.8 KiB
C++
105 lines
2.8 KiB
C++
/*
|
|
* Author: Andrej Prochazka
|
|
* Copyright(C) 2017
|
|
* This program is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU Affero General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Affero General Public License for more details.
|
|
|
|
* You should have received a copy of the GNU Affero General Public License
|
|
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#ifndef INTEGRATOR_H
|
|
#define INTEGRATOR_H
|
|
#include "catima/build_config.h"
|
|
#include "gsl/gsl_integration.h"
|
|
#include <functional>
|
|
#include <array>
|
|
#ifdef USE_THREADS
|
|
#include <mutex>
|
|
#endif
|
|
|
|
namespace catima{
|
|
|
|
/// helper class to integrate functions using the GSL library
|
|
class IntegratorGSL{
|
|
public:
|
|
IntegratorGSL(bool adapt=true);
|
|
~IntegratorGSL();
|
|
|
|
double integrate(std::function<double(double)> f, double min, double max, double precision=0.001);
|
|
|
|
private:
|
|
gsl_integration_workspace *w;
|
|
bool adaptive;
|
|
double error;
|
|
double result;
|
|
double min;
|
|
double max;
|
|
#ifdef USE_THREADS
|
|
std::mutex integration_mutex;
|
|
#endif
|
|
};
|
|
|
|
// built in integrator
|
|
template<int order>
|
|
class GaussLegendreIntegration{
|
|
public:
|
|
template<typename F>
|
|
double integrate(F f, double a, double b) const;
|
|
template<typename F>
|
|
double operator()(F f, double a, double b) const {return integrate(f, a, b);}
|
|
double get_w(int i) const {return w[i];}
|
|
double get_x(int i) const {return x[i];}
|
|
int n() const {return order;}
|
|
std::array<double,order> get_points(double a = -1.0, double b = 1.0)const;
|
|
|
|
public:
|
|
static std::array<double,order/2> w;
|
|
static std::array<double,order/2> x;
|
|
};
|
|
|
|
template<int order>
|
|
template<typename F>
|
|
double GaussLegendreIntegration<order>::integrate(F f, double a, double b) const{
|
|
double res=0.0;
|
|
double p = 0.5*(b-a);
|
|
double q = 0.5*(b+a);
|
|
|
|
for(int i=0;i<order/2;i++){
|
|
res += w[i] * (f(p*x[i] + q) + f(-p*x[i] + q));
|
|
}
|
|
return p*res;
|
|
}
|
|
|
|
template<int order>
|
|
std::array<double,order> GaussLegendreIntegration<order>::get_points(double a, double b)const{
|
|
std::array<double,order> points;
|
|
double p = 0.5*(b-a);
|
|
double q = 0.5*(b+a);
|
|
|
|
int num = (order/2);
|
|
for(int i=0;i< num;i++){
|
|
points[num-i-1] = -p*x[i] + q;
|
|
points[num+i] = p*x[i] + q;
|
|
}
|
|
return points;
|
|
}
|
|
|
|
#ifdef GSL_INTEGRATION
|
|
using integrator_type = IntegratorGSL;
|
|
#else
|
|
using integrator_type = GaussLegendreIntegration<8>;
|
|
#endif
|
|
|
|
extern integrator_type integrator;
|
|
extern IntegratorGSL integratorGSL;
|
|
}
|
|
|
|
#endif
|