Skip to content
GitLab
Commit a72fc0c3 authored by Joseph Boudreau's avatar Joseph Boudreau
Browse files

DefiniteIntegral now can handle improper integrals

parent 1dce6d26
Hide whitespace changes
Inline Side-by-side
GenericFunctions/GenericFunctions/DefiniteIntegral.hh
View file @ a72fc0c3
......@@ -2,9 +2,14 @@
// $Id: DefiniteIntegral.hh,v 1.2 2003/09/06 14:04:13 boudreau Exp $
//-------------------------------------------------------------//
// //
// This functional returns the definite integral of a function //
// This functional performs Romberg integration on a function //
// between lower bound and upper bound b. //
// //
// Two types: OPEN: use open quadrature formula //
// for improper integrals //
// CLOSED (default) use closed quadrature //
// formula. //
// //
//-------------------------------------------------------------//
#ifndef _DefiniteIntegral_h_
......@@ -21,8 +26,17 @@ namespace Genfun {
public:
// Type definition:
typedef enum {CLOSED, OPEN} Type;
// Constructor:
DefiniteIntegral(double a, double b);
DefiniteIntegral(double a, double b, Type=CLOSED);
// Copy Constructor:
DefiniteIntegral(const DefiniteIntegral &);
// Assignment Operator:
DefiniteIntegral & operator=(const DefiniteIntegral &) ;
// Destructor:
~DefiniteIntegral();
......@@ -33,23 +47,23 @@ namespace Genfun {
// Retrieve the number of function calls for the last operation:
unsigned int numFunctionCalls() const;
private:
// Algorithmic parameters:
// Desired precision (default 1.0E-06)
void setEpsilon(double eps);
// Trapezoid calculation:
double _trapzd( const AbsFunction & function, double a, double b, int j) const;
// Maximum number of iterations (default 20(closed) 14 (open))
void setMaxIter (unsigned int maxIter);
// Polynomial interpolation:
void _polint(double *xArray, double *yArray, double x, double & y, double & deltay) const;
// Minimum order:
void setMinOrder (unsigned int order);
double _a; // lower limit of integration
double _b; // upper limit of integration
static const int _K; // Order
static const int _KP; // Const dim of certain arrays.
// buffered value for _trapzd calculation:
mutable double _sTrap;
mutable unsigned int _nFunctionCalls;
private:
class Clockwork;
Clockwork *c;
};
} // namespace Genfun
#endif
GenericFunctions/src/DefiniteIntegral.cc
View file @ a72fc0c3
......@@ -2,111 +2,270 @@
// $Id: DefiniteIntegral.cc,v 1.6 2010/06/16 18:22:01 garren Exp $
#include <cmath>
#include <vector>
#include <stdexcept>
#include "CLHEP/GenericFunctions/DefiniteIntegral.hh"
#include "CLHEP/GenericFunctions/AbsFunction.hh"
namespace Genfun {
const int DefiniteIntegral::_K =5;
const int DefiniteIntegral::_KP=_K+1;
DefiniteIntegral::DefiniteIntegral(double a, double b):
_a(a), _b(b)
{}
DefiniteIntegral::~DefiniteIntegral() {
}
class DefiniteIntegral::Clockwork {
//
// This class has limited visibility its functions, data,
// and nested classes are all public:
//
public:
class QuadratureRule {
public:
// Constructorctor:
QuadratureRule() {};
// Destructor:
~QuadratureRule() {};
// Integrate at the j^th level of refinement:
virtual double integrate(const AbsFunction & function,
double a,
double b,
unsigned int j) const=0;
// Return the step multiplier:
virtual unsigned int stepMultiplier () const=0;
// Return the number of function calls:
virtual unsigned int numFunctionCalls() const =0;
};
class TrapezoidQuadratureRule:public QuadratureRule {
public:
// Constructor:
TrapezoidQuadratureRule():retVal(0),nFunctionCalls(0) {};
// Destructor:
~TrapezoidQuadratureRule() {};
// Integrate at the j^th level of refinement:
virtual double integrate(const AbsFunction & function,
double a,
double b,
unsigned int j) const;
// The step is doubled at each level of refinement:
virtual unsigned int stepMultiplier () const {return 2;}
#define FANCY
double DefiniteIntegral::operator [] (const AbsFunction & function) const {
_nFunctionCalls=0;
const int MAXITERATIONS=40; // Maximum number of iterations
const double EPS=1.0E-6; // Precision
#ifdef FANCY
// Returns number of function calls:
virtual unsigned int numFunctionCalls() const {return nFunctionCalls;};
private:
double s[MAXITERATIONS+2],h[MAXITERATIONS+2];
h[1]=1.0;
for (int j=1;j<=MAXITERATIONS;j++) {
s[j]=_trapzd(function, _a, _b, j);
if (j>=_K) {
double ss, dss;
_polint(h+j-_K,s+j-_K,0.0,ss, dss);
if (fabs(dss) <= EPS*fabs(ss)) return ss;
}
s[j+1]=s[j];
h[j+1]=0.25*h[j];
mutable double retVal;
mutable unsigned int nFunctionCalls;
};
class XtMidpointQuadratureRule:public QuadratureRule {
public:
// Constructor:
XtMidpointQuadratureRule():retVal(0),nFunctionCalls(0) {};
// Destructorctor:
~XtMidpointQuadratureRule() {};
// Integrate at the j^th level of refinement:
virtual double integrate(const AbsFunction & function,
double a,
double b,
unsigned int j) const;
// The step is tripled at each level of refinement:
virtual unsigned int stepMultiplier () const {return 3;}
// Returns number of function calls:
virtual unsigned int numFunctionCalls() const {return nFunctionCalls;};
private:
mutable double retVal;
mutable unsigned int nFunctionCalls;
};
double a; // lower limit of integration
double b; // upper limit of integration
DefiniteIntegral::Type type; // open or closed
mutable unsigned int nFunctionCalls; // bookkeeping
unsigned int MAXITER; // Max no of step doubling, tripling, etc.
double EPS; // Target precision
unsigned int K; // Minimum order =2*5=10
// Polynomial interpolation with Neville's method:
void polint(std::vector<double>::iterator xArray, std::vector<double>::iterator yArray, double x, double & y, double & deltay) const;
};
DefiniteIntegral::DefiniteIntegral(double a, double b, Type type):
c(new Clockwork()) {
c->a = a;
c->b = b;
c->type = type;
c->nFunctionCalls = 0;
c->MAXITER = type==OPEN ? 20 : 14;
c->EPS = 1.0E-6;
c->K = 5;
}
#else
double olds = -1E-30;
for (int j=1;j<MAXITERATIONS;j++) {
double s = _trapzd(function, _a,_b,j);
if (fabs(s-olds)<=EPS*fabs(olds)) return s;
olds=s;
DefiniteIntegral::~DefiniteIntegral() {
delete c;
}
#endif
std::cerr
<< "DefiniteIntegral: too many steps. No convergence"
<< std::endl;
return 0.0; // Never get here.
}
double DefiniteIntegral::_trapzd(const AbsFunction & function, double a, double b, int n) const {
int it, j;
if (n==1) {
_sTrap = 0.5*(b-a)*(function(a)+function(b));
_nFunctionCalls+=2;
DefiniteIntegral::DefiniteIntegral(const DefiniteIntegral & right) :c(new Clockwork(*right.c)) {
}
else {
for (it=1,j=1;j<n-1;j++) it <<=1;
double tnm=it;
double del = (b-a)/tnm;
double x=a+0.5*del;
double sum;
for (sum=0.0,j=1;j<=it;j++,x+=del) {
sum +=function(x);
_nFunctionCalls++;
DefiniteIntegral & DefiniteIntegral::operator = (const DefiniteIntegral & right) {
if (this!=&right) {
delete c;
c = new Clockwork(*right.c);
}
_sTrap = 0.5*(_sTrap+(b-a)*sum/tnm);
return *this;
}
void DefiniteIntegral::setEpsilon(double eps) {
c->EPS=eps;
}
return _sTrap;
}
void DefiniteIntegral::setMaxIter(unsigned int maxIter) {
c->MAXITER=maxIter;
}
void DefiniteIntegral::_polint(double *xArray, double *yArray, double x, double & y, double & deltay) const {
void DefiniteIntegral::setMinOrder(unsigned int minOrder) {
c->K=(minOrder+1)/2;
}
double dif = fabs(x-xArray[1]),dift;
double c[_KP],d[_KP];
int ns=1;
for (int i=1;i<=_K;i++) {
dift=fabs(x-xArray[i]);
if (dift<dif) {
ns=i;
dif=dift;
double DefiniteIntegral::operator [] (const AbsFunction & function) const {
const Clockwork::QuadratureRule * rule = c->type==OPEN ?
(Clockwork::QuadratureRule *) new Clockwork::XtMidpointQuadratureRule() :
(Clockwork::QuadratureRule *) new Clockwork::TrapezoidQuadratureRule();
double xMult=rule->stepMultiplier();
c->nFunctionCalls=0;
std::vector<double> s(c->MAXITER+2),h(c->MAXITER+2);
h[1]=1.0;
for (unsigned int j=1;j<=c->MAXITER;j++) {
s[j]=rule->integrate(function, c->a, c->b, j);
c->nFunctionCalls=rule->numFunctionCalls();
if (j>=c->K) {
double ss, dss;
c->polint(h.begin()+j-c->K,s.begin()+j-c->K,0.0,ss, dss);
if (fabs(dss) <= c->EPS*fabs(ss)) {
delete rule;
return ss;
}
}
s[j+1]=s[j];
h[j+1]=h[j]/xMult/xMult;
}
c[i]=d[i]=yArray[i];
delete rule;
throw std::runtime_error("DefiniteIntegral: too many steps. No convergence");
return 0.0; // Never get here.
}
y = yArray[ns--];
for (int m=1;m<_K;m++) {
for (int i=1;i<=_K-m;i++) {
double ho = xArray[i]-x;
double hp= xArray[i+m]-x;
double w=c[i+1]-d[i];
double den=ho-hp;
if (den==0)
std::cerr
<< "Error in polynomial extrapolation"
<< std::endl;
den=w/den;
d[i]=hp*den;
c[i]=ho*den;
void DefiniteIntegral::Clockwork::polint(std::vector<double>::iterator xArray, std::vector<double>::iterator yArray, double x, double & y, double & deltay) const {
double dif = fabs(x-xArray[1]),dift;
std::vector<double> c(K+1),d(K+1);
int ns=1;
for (unsigned int i=1;i<=K;i++) {
dift=fabs(x-xArray[i]);
if (dift<dif) {
ns=i;
dif=dift;
}
c[i]=d[i]=yArray[i];
}
y = yArray[ns--];
for (unsigned int m=1;m<K;m++) {
for (unsigned int i=1;i<=K-m;i++) {
double ho = xArray[i]-x;
double hp= xArray[i+m]-x;
double w=c[i+1]-d[i];
double den=ho-hp;
if (den==0)
std::cerr
<< "Error in polynomial extrapolation"
<< std::endl;
den=w/den;
d[i]=hp*den;
c[i]=ho*den;
}
deltay = 2*ns < (K-m) ? c[ns+1]: d[ns--];
y += deltay;
}
deltay = 2*ns < (_K-m) ? c[ns+1]: d[ns--];
y += deltay;
}
}
unsigned int DefiniteIntegral::numFunctionCalls() const {
return _nFunctionCalls;
return c->nFunctionCalls;
}
// Quadrature rules:
double DefiniteIntegral::Clockwork::TrapezoidQuadratureRule::integrate(const AbsFunction & function, double a, double b, unsigned int n) const {
unsigned int it, j;
if (n==1) {
retVal = 0.5*(b-a)*(function(a)+function(b));
nFunctionCalls+=2;
}
else {
for (it=1,j=1;j<n-1;j++) it <<=1;
double tnm=it;
double del = (b-a)/tnm;
double x=a+0.5*del;
double sum;
for (sum=0.0,j=1;j<=it;j++,x+=del) {
sum +=function(x);
nFunctionCalls++;
}
retVal = 0.5*(retVal+(b-a)*sum/tnm);
}
return retVal;
}
// Quadrature rules:
double DefiniteIntegral::Clockwork::XtMidpointQuadratureRule::integrate(const AbsFunction & function, double a, double b, unsigned int n) const {
unsigned int it, j;
if (n==1) {
retVal = (b-a)*(function((a+b)/2.0));
nFunctionCalls+=1;
}
else {
for (it=1,j=1;j<n-1;j++) it *=3;
double tnm=it;
double del = (b-a)/(3.0*tnm);
double ddel = del+del;
double x=a+0.5*del;
double sum=0;
for (j=1;j<=it;j++) {
sum +=function(x);
x+=ddel;
sum +=function(x);
x+=del;
nFunctionCalls+=2;
}
retVal = (retVal+(b-a)*sum/tnm)/3.0;
}
return retVal;
}
} // namespace Genfun
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