#include #include #include #include #include "JTools/JQuadrature.hh" #include "JMath/JPolynome.hh" #include "JMath/JGauss.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace GAUSS_LEGENDRE { using namespace JPP; const int numberOfPoints = 10; const double epsilon = 1.0e-12; const double xmin = -5.0; const double xmax = +3.0; /** * Function. */ const JPolynome f1(1.0, 2.0, 3.0); /** * Integral. */ const JPolynome F1 = f1.getIntegral(); } namespace GAUSS_LAGUERRE { using namespace JPP; const int numberOfPoints = 10; const double a = 1.0; const double epsilon = 1.0e-12; const double xmin = 0.0; const double xmax = +1.0e10; // approximation of infinity /** * Function. */ inline double f1(const double x) { return sin(x); } /** * Integral. */ inline double F1(const double x) { // integral of sin(x) * x^1 * exp(-x) return -0.5 * exp(-x) * (x*sin(x) + (x + 1.0)*cos(x)); } } namespace GAUSS_HERMITE { using namespace JPP; const int numberOfPoints = 10; const double epsilon = 1.0e-12; const double sigma[] = { 4.0, 0.5 }; const double xmin = -3.0 * sigma[0]; const double xmax = +3.0 * sigma[0]; /** * Function. */ const JGauss f1(0.0, sigma[0]); /** * Integral. */ const JGauss F1(0.0, hypot(sigma[0], sigma[1])); } /** * \file * * Example program to test JTOOLS::JQuadrature. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; double precision; int debug; try { JParser<> zap("Example program to test JTOOLS::JQuadrature."); zap['e'] = make_field(precision) = 1.0e-4; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } { using namespace GAUSS_LEGENDRE; double V = F1(xmax) - F1(xmin); JGaussLegendre engine(numberOfPoints, epsilon); const double x0 = 0.5 * (xmax + xmin); const double x1 = 0.5 * (xmax - xmin); double W = 0.0; for (JGaussLegendre::const_iterator i = engine.begin(); i != engine.end(); ++i) { const double x = x0 + i->getX() * x1; W += f1(x) * i->getY() * x1; } DEBUG("integral analytical " << FIXED(9,5) << V << endl); DEBUG("integral numerical " << FIXED(9,5) << W << endl); ASSERT(fabs(V - W) <= precision, "Test Gauss-Legendre integration."); } { using namespace GAUSS_LAGUERRE; const double V = F1(xmax) - F1(xmin); JGaussLaguerre engine(numberOfPoints, a, epsilon); double W = 0.0; for (JGaussLaguerre::const_iterator i = engine.begin(); i != engine.end(); ++i) { const double x = i->getX(); W += f1(x) * i->getY(); } DEBUG("integral analytical " << FIXED(9,5) << V << endl); DEBUG("integral numerical " << FIXED(9,5) << W << endl); ASSERT(fabs(V - W) <= precision, "Test Gauss-Laguerre integration."); } { using namespace GAUSS_HERMITE; JGaussHermite engine(numberOfPoints, epsilon); for (double x = xmin; x <= xmax; x += 0.05*(xmax - xmin)) { const double V = F1(x); double W = 0.0; for (JGaussHermite::const_iterator i = engine.begin(); i != engine.end(); ++i) { const double u = i->getX() * sigma[1] * sqrt(2.0); const double v = i->getY() / sqrt(PI); W += f1(x + u) * v; } DEBUG(FIXED(7,3) << x << ' ' << FIXED(9,5) << V << ' ' << FIXED(9,5) << W << endl); ASSERT(fabs(V - W) <= precision * V, "Test Gauss-Hermite integration."); } } return 0; }