JQuadrature.cc File Reference

Example program to test JTOOLS::JQuadrature. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "JTools/JQuadrature.hh"
#include "JMath/JPolynome.hh"
#include "JMath/JGauss.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Namespaces
	GAUSS_LEGENDRE
	GAUSS_LAGUERRE
	GAUSS_HERMITE

Functions
const JPolynome	GAUSS_LEGENDRE::f1 (1.0, 2.0, 3.0)
	Function. More...
double	GAUSS_LAGUERRE::f1 (const double x)
	Function. More...
double	GAUSS_LAGUERRE::F1 (const double x)
	Integral. More...
const JGauss	GAUSS_HERMITE::f1 (0.0, sigma[0])
	Function. More...
const JGauss	GAUSS_HERMITE::F1 (0.0, hypot(sigma[0], sigma[1]))
	Integral. More...
int	main (int argc, char **argv)

Variables
const int	GAUSS_LEGENDRE::numberOfPoints = 10
const double	GAUSS_LEGENDRE::epsilon = 1.0e-12
const double	GAUSS_LEGENDRE::xmin = -5.0
const double	GAUSS_LEGENDRE::xmax = +3.0
const JPolynome	GAUSS_LEGENDRE::F1 = f1.getIntegral()
	Integral. More...
const int	GAUSS_LAGUERRE::numberOfPoints = 10
const double	GAUSS_LAGUERRE::a = 1.0
const double	GAUSS_LAGUERRE::epsilon = 1.0e-12
const double	GAUSS_LAGUERRE::xmin = 0.0
const double	GAUSS_LAGUERRE::xmax = +1.0e10
const int	GAUSS_HERMITE::numberOfPoints = 10
const double	GAUSS_HERMITE::epsilon = 1.0e-12
const double	GAUSS_HERMITE::sigma [] = { 4.0, 0.5 }
const double	GAUSS_HERMITE::xmin = -3.0 * sigma[0]
const double	GAUSS_HERMITE::xmax = +3.0 * sigma[0]

Detailed Description

Example program to test JTOOLS::JQuadrature.

Author

mdejong

Definition in file JQuadrature.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 97 of file JQuadrature.cc.

{
  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;
}