Commit cc22ab32 authored by Pawel Sznajder's avatar Pawel Sznajder
Browse files

Add new integration routine TRAPEIZODALLOG with trapeizodal integration

rule probing function in equidistant steps between ln(a) and ln(b) 
parent b1487da8
......@@ -26,7 +26,7 @@ public:
* Type of integration.
*/
enum Type {
UNDEFINED = 0, TRAPEZOIDAL, GaussLegendre, DEXP, GK21_ADAPTIVE
UNDEFINED = 0, TRAPEZOIDAL, TRAPEZOIDALLOG, GaussLegendre, DEXP, GK21_ADAPTIVE
};
/**
......
......@@ -9,7 +9,6 @@
* @version 1.0
*/
//#include <string>
#include "QuadratureIntegrator1D.h"
class FunctionType1D;
......@@ -25,12 +24,15 @@ namespace NumA {
*/
class TrapezoidalIntegrator1D: public QuadratureIntegrator1D {
public:
/**
* Constructor.
* @param N Order of the quadrature (number of nodes).
*/
TrapezoidalIntegrator1D(unsigned int N = 50);
/**
* Default destructor.
*/
......@@ -40,18 +42,19 @@ public:
virtual void setN(unsigned int n);
// Not needed anymore, done in parent class.
// virtual double integrate(FunctionType1D* pFunction, double a, double b,
// std::vector<double>& parameters);
protected:
/**
* Copy constructor.
* Called by clone().
* @param other TrapezoidalIntegrator1D object to copy.
*/
TrapezoidalIntegrator1D(const TrapezoidalIntegrator1D &other);
void setNodesAndWeights(); ///< Computes the nodes and weights of the Gauss-Legendre quadrature. Called when N is set.
/**
* Computes the nodes and weights of the Gauss-Legendre quadrature. Called when N is set.
*/
void setNodesAndWeights();
};
} /* namespace NumA */
......
#ifndef TRAPEZOIDAL_LOG_INTEGRATOR_1D_H
#define TRAPEZOIDAL_LOG_INTEGRATOR_1D_H
/**
* @file TrapezoidalLogIntegrator1D.h
* @author Pawel Sznajder (NCBJ, Warsaw)
* @date 12 July 2017
* @version 1.0
*/
#include <ElementaryUtils/parameters/Parameters.h>
#include <include/NumA/integration/one_dimension/QuadratureIntegrator1D.h>
#include <string>
#include <vector>
class FunctionType1D;
namespace NumA {
/**
* @class TrapezoidalLogIntegrator1D
*
* @brief Trapezoidal integration with logarithmic step.
*
* See Integrator1D documentation for an example.
*/
class TrapezoidalLogIntegrator1D: public QuadratureIntegrator1D {
public:
/**
* Constructor.
* @param N Order of the quadrature (number of nodes).
*/
TrapezoidalLogIntegrator1D();
/**
* Default destructor.
*/
virtual ~TrapezoidalLogIntegrator1D();
virtual TrapezoidalLogIntegrator1D* clone() const;
virtual void configure(const ElemUtils::Parameters &parameters);
virtual void setN(unsigned int n);
virtual double integrate(FunctionType1D* pFunction, double a, double b,
std::vector<double> &parameters);
protected:
/**
* Copy constructor.
* Called by clone().
* @param other TrapezoidalLogIntegrator1D object to copy.
*/
TrapezoidalLogIntegrator1D(const TrapezoidalLogIntegrator1D &other);
};
} /* namespace NumA */
#endif /* TRAPEZOIDAL_LOG_INTEGRATOR_1D_H */
......@@ -10,6 +10,7 @@
#include "../../../../include/NumA/integration/one_dimension/GaussLegendreIntegrator1D.h"
#include "../../../../include/NumA/integration/one_dimension/IntegratorType1D.h"
#include "../../../../include/NumA/integration/one_dimension/TrapezoidalIntegrator1D.h"
#include "../../../../include/NumA/integration/one_dimension/TrapezoidalLogIntegrator1D.h"
#include "../../../../include/NumA/utils/Errors.h"
#include "../../../../include/NumA/utils/Tolerances.h"
......@@ -37,6 +38,9 @@ Integrator1D* Integrator1D::newIntegrator(
case IntegratorType1D::TRAPEZOIDAL:
return new TrapezoidalIntegrator1D();
break;
case IntegratorType1D::TRAPEZOIDALLOG:
return new TrapezoidalLogIntegrator1D();
break;
case IntegratorType1D::GaussLegendre:
return new GaussLegendreIntegrator1D();
break;
......
......@@ -22,6 +22,9 @@ std::string IntegratorType1D::toString() {
case TRAPEZOIDAL:
return "TRAPEZOIDAL";
break;
case TRAPEZOIDALLOG:
return "TRAPEZOIDALLOG";
break;
case GaussLegendre:
return "Gauss-Legendre";
break;
......@@ -41,6 +44,9 @@ std::string IntegratorType1D::getShortName() {
case TRAPEZOIDAL:
return "TRAPEZOIDAL";
break;
case TRAPEZOIDALLOG:
return "TRAPEZOIDALLOG";
break;
case GaussLegendre:
return "GL";
break;
......
......@@ -2,8 +2,6 @@
#include <vector>
//#include "../../../../include/NumA/utils/Parameters.h"
namespace NumA {
TrapezoidalIntegrator1D::TrapezoidalIntegrator1D(unsigned int N) :
......@@ -46,26 +44,4 @@ void TrapezoidalIntegrator1D::setNodesAndWeights() {
}
}
// Not needed anymore, done in parent class.
//double TrapezoidalIntegrator1D::integrate(FunctionType1D* pFunction, double a,
// double b, std::vector<double>& parameters) {
// double result = 0., x = 0.;
// // Defining integration step
// double step = (b - a) / (m_N - 1);
//
// // Adding the border values before integrating through the loop
// result = 0.5 * (*pFunction)(a, parameters);
// result += 0.5 * (*pFunction)(b, parameters);
//
// // Entering the loop
// for (size_t i = 1; i < m_N - 1; i++) {
// x = a + i * step;
// result += (*pFunction)(x, parameters);
// }
// // Multiplying by the step at the end to save computing time
// result = result * step;
//
// return result;
//}
} /* namespace NumA */
#include <ElementaryUtils/logger/CustomException.h>
#include <ElementaryUtils/string_utils/Formatter.h>
#include <include/NumA/functor/one_dimension/FunctionType1D.h>
#include <include/NumA/integration/one_dimension/TrapezoidalLogIntegrator1D.h>
#include <stddef.h>
#include <cmath>
namespace NumA {
TrapezoidalLogIntegrator1D::TrapezoidalLogIntegrator1D() :
QuadratureIntegrator1D() {
}
TrapezoidalLogIntegrator1D::TrapezoidalLogIntegrator1D(
const TrapezoidalLogIntegrator1D &other) :
QuadratureIntegrator1D(other) {
}
TrapezoidalLogIntegrator1D::~TrapezoidalLogIntegrator1D() {
}
TrapezoidalLogIntegrator1D* TrapezoidalLogIntegrator1D::clone() const {
return new TrapezoidalLogIntegrator1D(*this);
}
double TrapezoidalLogIntegrator1D::integrate(FunctionType1D* pFunction,
double a, double b, std::vector<double> &parameters) {
//check if valid
if (a <= 0. || b <= 0.) {
throw ElemUtils::CustomException("TrapezoidalLogIntegrator1D", __func__,
ElemUtils::Formatter()
<< "This integration routine can be used only for > 0 range, here: "
<< a << ", " << b);
}
//result
double result = 0.;
//log range and step
double logA = log(a);
double logB = log(b);
double logStep = (logB - logA) / double(m_N - 1);
double argumentA = a;
double argumentB;
double valueA = ((*pFunction)(a, parameters));
double valueB;
// Entering the loop
for (size_t i = 1; i < m_N; i++) {
//actual
argumentB = exp(logA + i * logStep);
//calculate
valueB = ((*pFunction)(argumentB, parameters));
//add
result += (argumentB - argumentA) * (valueA + valueB);
//save previous
argumentA = argumentB;
valueA = valueB;
}
//return
return 0.5 * result;
}
void TrapezoidalLogIntegrator1D::configure(
const ElemUtils::Parameters &parameters) {
QuadratureIntegrator1D::configure(parameters);
}
void TrapezoidalLogIntegrator1D::setN(unsigned int n) {
m_N = n;
}
} /* namespace NumA */
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment