Commit 14210754 authored by Daniele Binosi's avatar Daniele Binosi
Browse files

ref#16

New classes that generates quadratures for numerical integration.

Chebyshev A: Chebyshev quadrature for computing the integral \int_{-1}^1 dx f(x)sqrt(1-x^2);
Chebyshev B: Chebyshev quadrature for computing the integral \int_{-1}^1 dx f(x)/sqrt(1-x^2);
GLSeSt: Different (more precise?) method for computing GL weights. 
parent 0f9da999
#ifndef CHEBYSHEVA_INTEGRATOR_1D_H
#define CHEBYSHEVA_INTEGRATOR_1D_H
/**
* @file ChebyshevAIntegrator1D.h
* @author Daniele BINOSI (ECT*)
* @date 27 April 2016
* @version 1.0
*/
//#include <string>
#include "QuadratureIntegrator1D.h"
class FunctionType1D;
namespace NumA {
/**
* @class ChebyshevAIntegrator1D
*
* @brief Chebyshev quadrature for computing the integral \int_{-1}^1 dx f(x)sqrt(1-x^2).
*/
class ChebyshevAIntegrator1D: public QuadratureIntegrator1D {
public:
ChebyshevAIntegrator1D(unsigned int N = 48);
virtual ~ChebyshevAIntegrator1D();
virtual ChebyshevAIntegrator1D* clone() const;
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:
ChebyshevAIntegrator1D(const ChebyshevAIntegrator1D &other);
void setNodesAndWeights();
};
} /* namespace NumA */
#endif /* CHEBYSHEVA_INTEGRATOR_1D_H */
#ifndef CHEBYSHEVB_INTEGRATOR_1D_H
#define CHEBYSHEVB_INTEGRATOR_1D_H
/**
* @file ChebyshevBIntegrator1D.h
* @author Daniele BINOSI (ECT*)
* @date 27 April 2016
* @version 1.0
*/
//#include <string>
#include "QuadratureIntegrator1D.h"
class FunctionType1D;
namespace NumA {
/**
* @class ChebyshevBIntegrator1D
*
* @brief Chebyshev quadrature for computing the integral \int_{-1}^1 dx f(x)/sqrt(1-x^2).
*/
class ChebyshevBIntegrator1D: public QuadratureIntegrator1D {
public:
ChebyshevBIntegrator1D(unsigned int N = 48);
virtual ~ChebyshevBIntegrator1D();
virtual ChebyshevBIntegrator1D* clone() const;
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:
ChebyshevBIntegrator1D(const ChebyshevBIntegrator1D &other);
void setNodesAndWeights();
};
} /* namespace NumA */
#endif /* CHEBYSHEVB_INTEGRATOR_1D_H */
#ifndef GAUSSLEGENDRESEST_INTEGRATOR_1D_H_
#define GAUSSLEGENDRESEST_INTEGRATOR_1D_H_
/**
* @file GaussLegendreSeStIntegrator1D.h
* @author Daniele BINOSI (ECT*)
* @date 6 May 2016
* @version 1.0
*/
//#include <string>
#include "QuadratureIntegrator1D.h"
namespace NumA {
/**
* @class GaussLegendreSestIntegrator1D
*
* @brief Gauss-Legendre quadrature.
*/
class GaussLegendreSeStIntegrator1D: public QuadratureIntegrator1D {
public:
GaussLegendreSeStIntegrator1D(unsigned int N = 50); // TODO how many steps by default ?
virtual ~GaussLegendreSeStIntegrator1D();
virtual void setN(unsigned int n);
virtual GaussLegendreSeStIntegrator1D* clone() const;
protected:
GaussLegendreSeStIntegrator1D(const GaussLegendreSeStIntegrator1D &other);
void setNodesAndWeights();
};
} /* namespace NumA */
#endif /* GAUSSLEGENDRE_INTEGRATOR_1D_H_ */
#include "../../../../include/NumA/integration/one_dimension/ChebyshevAIntegrator1D.h"
#include <stdexcept>
#include <vector>
#include <math.h>
namespace NumA {
ChebyshevAIntegrator1D::ChebyshevAIntegrator1D(unsigned int N) :
QuadratureIntegrator1D(N) {
setNodesAndWeights();
}
ChebyshevAIntegrator1D::ChebyshevAIntegrator1D(
const ChebyshevAIntegrator1D &other) :
QuadratureIntegrator1D(other) {
}
ChebyshevAIntegrator1D::~ChebyshevAIntegrator1D() {
}
ChebyshevAIntegrator1D* ChebyshevAIntegrator1D::clone() const {
return new ChebyshevAIntegrator1D(*this);
}
void ChebyshevAIntegrator1D::setN(unsigned int n) {
if (n != m_N) {
m_N = n;
setNodesAndWeights();
}
}
void ChebyshevAIntegrator1D::setNodesAndWeights() {
double phi = 0.;
//clear nodes and weights
m_nodes.clear();
m_weights.clear();
// Check number of steps
if (m_N < 1) {
throw std::runtime_error(
"[ChebyshevAIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 1");
} else {
m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.);
for (unsigned int i = 0; i <m_N; i++) {
phi = M_PI*double(m_N-i)/double(m_N+1);
m_nodes.at(i) = cos(phi);
m_weights.at(i) = M_PI/double(m_N+1)*(sin(phi)*sin(phi));
}
}
}
} /* namespace NumA */
#include "../../../../include/NumA/integration/one_dimension/ChebyshevBIntegrator1D.h"
#include <stdexcept>
#include <vector>
#include <math.h>
namespace NumA {
ChebyshevBIntegrator1D::ChebyshevBIntegrator1D(unsigned int N) :
QuadratureIntegrator1D(N) {
setNodesAndWeights();
}
ChebyshevBIntegrator1D::ChebyshevBIntegrator1D(
const ChebyshevBIntegrator1D &other) :
QuadratureIntegrator1D(other) {
}
ChebyshevBIntegrator1D::~ChebyshevBIntegrator1D() {
}
ChebyshevBIntegrator1D* ChebyshevBIntegrator1D::clone() const {
return new ChebyshevBIntegrator1D(*this);
}
void ChebyshevBIntegrator1D::setN(unsigned int n) {
if (n != m_N) {
m_N = n;
setNodesAndWeights();
}
}
void ChebyshevBIntegrator1D::setNodesAndWeights() {
double phi = 0.;
//clear nodes and weights
m_nodes.clear();
m_weights.clear();
// Check number of steps
if (m_N < 1) {
throw std::runtime_error(
"[ChebyshevBIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 1");
} else {
m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.);
for (unsigned int i = 0; i < m_N; i++) {
phi = M_PI*double(2*m_N-1-2*i)/double(2*m_N);
m_nodes.at(i) = cos(phi);
m_weights.at(i) = M_PI/double(m_N);
}
}
}
} /* namespace NumA */
#include "../../../../include/NumA/integration/one_dimension/GaussLegendreSeStIntegrator1D.h"
#include <algorithm>
#include <cmath>
//#include <functional>
//#include <iostream>
#include <limits>
//#include <numeric>
//#include <sstream>
#include <stdexcept>
#include <vector>
//#include "../../../../include/NumA/utils/GenericType.h"
//#include "../../../../include/NumA/utils/Parameters.h"
namespace NumA {
GaussLegendreSeStIntegrator1D::GaussLegendreSeStIntegrator1D(unsigned int N) :
QuadratureIntegrator1D(N) {
setNodesAndWeights();
}
GaussLegendreSeStIntegrator1D::GaussLegendreSeStIntegrator1D(
const GaussLegendreSeStIntegrator1D &other) :
QuadratureIntegrator1D(other) {
}
GaussLegendreSeStIntegrator1D::~GaussLegendreSeStIntegrator1D() {
// TODO Auto-generated destructor stub
}
GaussLegendreSeStIntegrator1D* GaussLegendreSeStIntegrator1D::clone() const {
return new GaussLegendreSeStIntegrator1D(*this);
}
void GaussLegendreSeStIntegrator1D::setN(unsigned int n) {
if (n != m_N) {
m_N = n;
setNodesAndWeights();
}
}
void GaussLegendreSeStIntegrator1D::setNodesAndWeights() {
//clear nodes and weights
m_nodes.clear();
m_weights.clear();
// Check number of steps
if (m_N <= 0) {
throw std::runtime_error(
"[GaussLegendreSeStIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 0");
} else {
m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.);
std::vector< double > c(m_N,0.);
double cc = 1.;
double r = 0.;
double xtemp = 0.;
unsigned int const maxit = 10;
double p0, dp0, p1, dp1, p2, dp2, d;
for (unsigned int i = 0; i <m_N; i++) {
c[i] = double(i*i)/double((2*i+1)*(2*i-1));
}
for (unsigned int i = 1; i <m_N; i++) {
cc *= c[i];
}
cc = 2*cc;
for (unsigned int i = 0; i <m_N; i++) {
if (i == 0) {
r = 2.78/(4.0+double(m_N*m_N));
xtemp = 1.0-r;
} else if (i == 1) {
r = 1.0+0.06*((double) m_N-8)/((double) m_N);
xtemp = xtemp-4.1*r*(1.0-xtemp);
} else if (i == 2) {
r = 1.0+0.22*((double) m_N-8)/((double) m_N);
xtemp = xtemp-1.67*r*(m_nodes[0]-xtemp);
} else if (i < m_N-2) {
xtemp = 3.0*m_nodes[i-1]-3.0*m_nodes[i-2]+m_nodes[i-3];
} else if (i == m_N-2) {
r = 1.0/(1.0+0.639*((double) m_N-4)/(1.0+0.71*((double) m_N-4)));
xtemp = xtemp+r*(xtemp-m_nodes[i-2])/0.766;
} else if (i == m_N-1) {
r = 1.0/(1.0+0.22*((double) m_N-8)/((double) m_N));
xtemp = xtemp+r*(xtemp-m_nodes[i-2])/1.67;
}
for (unsigned int j = 1; j <= maxit; j++) {
p1 = 1.;
dp1 = 0.;
p2 = xtemp;
dp2 = 1.;
for (unsigned int k = 1; k < m_N; k++) {
p0 = p1;
dp0 = dp1;
p1 = p2;
dp1 = dp2;
p2 = xtemp*p1-c[k]*p0;
dp2 = xtemp*dp1+p1-c[k]*dp0;
}
d = p2/dp2;
xtemp = xtemp-d;
if (std::abs(d) < std::numeric_limits<double>::epsilon()*(std::abs(d)+1.)) break;
}
m_nodes[i] = xtemp;
m_weights[i] = cc/dp2/p1;
}
std::reverse(m_nodes.begin(), m_nodes.end());
std::reverse(m_weights.begin(), m_weights.end());
}
}
} /* namespace NumA */
Supports Markdown
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