CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

DefiniteIntegral.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: DefiniteIntegral.cc,v 1.6 2010/06/16 18:22:01 garren Exp $
3 
4 #include <cmath>
5 #include <vector>
6 #include <stdexcept>
9 
10 
11 namespace Genfun {
12 
13 
15 
16  //
17  // This class has limited visibility its functions, data,
18  // and nested classes are all public:
19  //
20 
21  public:
22 
24 
25  public:
26 
27  // Constructorctor:
29 
30  // Destructor:
31  virtual ~QuadratureRule() {};
32 
33  // Integrate at the j^th level of refinement:
34  virtual double integrate(const AbsFunction & function,
35  double a,
36  double b,
37  unsigned int j) const=0;
38 
39  // Return the step multiplier:
40  virtual unsigned int stepMultiplier () const=0;
41 
42  // Return the number of function calls:
43  virtual unsigned int numFunctionCalls() const =0;
44 
45  };
46 
48 
49  public:
50 
51  // Constructor:
52  TrapezoidQuadratureRule():retVal(0),nFunctionCalls(0) {};
53 
54  // Destructor:
56 
57  // Integrate at the j^th level of refinement:
58  virtual double integrate(const AbsFunction & function,
59  double a,
60  double b,
61  unsigned int j) const;
62 
63  // The step is doubled at each level of refinement:
64  virtual unsigned int stepMultiplier () const {return 2;}
65 
66  // Returns number of function calls:
67  virtual unsigned int numFunctionCalls() const {return nFunctionCalls;};
68 
69  private:
70 
71  mutable double retVal;
72  mutable unsigned int nFunctionCalls;
73 
74  };
75 
77 
78  public:
79 
80  // Constructor:
81  XtMidpointQuadratureRule():retVal(0),nFunctionCalls(0) {};
82 
83  // Destructorctor:
85 
86  // Integrate at the j^th level of refinement:
87  virtual double integrate(const AbsFunction & function,
88  double a,
89  double b,
90  unsigned int j) const;
91 
92  // The step is tripled at each level of refinement:
93  virtual unsigned int stepMultiplier () const {return 3;}
94 
95  // Returns number of function calls:
96  virtual unsigned int numFunctionCalls() const {return nFunctionCalls;};
97 
98  private:
99 
100  mutable double retVal;
101  mutable unsigned int nFunctionCalls;
102  };
103 
104  double a; // lower limit of integration
105  double b; // upper limit of integration
106  DefiniteIntegral::Type type; // open or closed
107  mutable unsigned int nFunctionCalls; // bookkeeping
108  unsigned int MAXITER; // Max no of step doubling, tripling, etc.
109  double EPS; // Target precision
110  unsigned int K; // Minimum order =2*5=10
111 
112  // Polynomial interpolation with Neville's method:
113  void polint(std::vector<double>::iterator xArray, std::vector<double>::iterator yArray, double x, double & y, double & deltay) const;
114 
115  };
116 
117 
119  c(new Clockwork()) {
120  c->a = a;
121  c->b = b;
122  c->type = type;
123  c->nFunctionCalls = 0;
124  c->MAXITER = type==OPEN ? 20 : 14;
125  c->EPS = 1.0E-6;
126  c->K = 5;
127  }
128 
130  delete c;
131  }
132 
134  :AbsFunctional(right), c(new Clockwork(*right.c)) {
135  }
136 
138  if (this!=&right) {
139  delete c;
140  c = new Clockwork(*right.c);
141  }
142  return *this;
143  }
144 
145  void DefiniteIntegral::setEpsilon(double eps) {
146  c->EPS=eps;
147  }
148 
149  void DefiniteIntegral::setMaxIter(unsigned int maxIter) {
150  c->MAXITER=maxIter;
151  }
152 
153  void DefiniteIntegral::setMinOrder(unsigned int minOrder) {
154  c->K=(minOrder+1)/2;
155  }
156 
157  double DefiniteIntegral::operator [] (const AbsFunction & function) const {
158 
159  const Clockwork::QuadratureRule * rule = c->type==OPEN ?
162  double xMult=rule->stepMultiplier();
163 
164  c->nFunctionCalls=0;
165  std::vector<double> s(c->MAXITER+2),h(c->MAXITER+2);
166  h[1]=1.0;
167  for (unsigned int j=1;j<=c->MAXITER;j++) {
168  s[j]=rule->integrate(function, c->a, c->b, j);
169  c->nFunctionCalls=rule->numFunctionCalls();
170  if (j>=c->K) {
171  double ss, dss;
172  c->polint(h.begin()+j-c->K,s.begin()+j-c->K,0.0,ss, dss);
173  if (fabs(dss) <= c->EPS*fabs(ss)) {
174  delete rule;
175  return ss;
176  }
177  }
178  s[j+1]=s[j];
179  h[j+1]=h[j]/xMult/xMult;
180  }
181  delete rule;
182  throw std::runtime_error("DefiniteIntegral: too many steps. No convergence");
183  return 0.0; // Never get here.
184  }
185 
186  void DefiniteIntegral::Clockwork::polint(std::vector<double>::iterator xArray, std::vector<double>::iterator yArray, double x, double & y, double & deltay) const {
187  double dif = fabs(x-xArray[1]),dift;
188  std::vector<double> c(K+1),d(K+1);
189  unsigned int ns=1;
190  for (unsigned int i=1;i<=K;i++) {
191  dift=fabs(x-xArray[i]);
192  if (dift<dif) {
193  ns=i;
194  dif=dift;
195  }
196  c[i]=d[i]=yArray[i];
197  }
198  y = yArray[ns--];
199  for (unsigned int m=1;m<K;m++) {
200  for (unsigned int i=1;i<=K-m;i++) {
201  double ho = xArray[i]-x;
202  double hp= xArray[i+m]-x;
203  double w=c[i+1]-d[i];
204  double den=ho-hp;
205  if (den==0)
206  std::cerr
207  << "Error in polynomial extrapolation"
208  << std::endl;
209  den=w/den;
210  d[i]=hp*den;
211  c[i]=ho*den;
212  }
213  deltay = 2*ns < (K-m) ? c[ns+1]: d[ns--];
214  y += deltay;
215  }
216  }
217 
218  unsigned int DefiniteIntegral::numFunctionCalls() const {
219  return c->nFunctionCalls;
220  }
221 
222  // Quadrature rules:
223  double DefiniteIntegral::Clockwork::TrapezoidQuadratureRule::integrate(const AbsFunction & function, double a, double b, unsigned int n) const {
224  unsigned int it, j;
225  if (n==1) {
226  retVal = 0.5*(b-a)*(function(a)+function(b));
227  nFunctionCalls+=2;
228  }
229  else {
230  for (it=1,j=1;j<n-1;j++) it <<=1;
231  double tnm=it;
232  double del = (b-a)/tnm;
233  double x=a+0.5*del;
234  double sum;
235  for (sum=0.0,j=1;j<=it;j++,x+=del) {
236  sum +=function(x);
237  nFunctionCalls++;
238  }
239  retVal = 0.5*(retVal+(b-a)*sum/tnm);
240  }
241  return retVal;
242  }
243 
244  // Quadrature rules:
245  double DefiniteIntegral::Clockwork::XtMidpointQuadratureRule::integrate(const AbsFunction & function, double a, double b, unsigned int n) const {
246  unsigned int it, j;
247  if (n==1) {
248  retVal = (b-a)*(function((a+b)/2.0));
249  nFunctionCalls+=1;
250  }
251  else {
252  for (it=1,j=1;j<n-1;j++) it *=3;
253  double tnm=it;
254  double del = (b-a)/(3.0*tnm);
255  double ddel = del+del;
256  double x=a+0.5*del;
257  double sum=0;
258  for (j=1;j<=it;j++) {
259  sum +=function(x);
260  x+=ddel;
261  sum +=function(x);
262  x+=del;
263  nFunctionCalls+=2;
264  }
265  retVal = (retVal+(b-a)*sum/tnm)/3.0;
266  }
267  return retVal;
268  }
269 
270 
271 
272 } // namespace Genfun
void setMaxIter(unsigned int maxIter)
DefiniteIntegral & operator=(const DefiniteIntegral &)
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const =0
void setMinOrder(unsigned int order)
virtual unsigned int numFunctionCalls() const =0
virtual double operator[](const AbsFunction &function) const
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const
virtual unsigned int stepMultiplier() const =0
DefiniteIntegral(double a, double b, Type=CLOSED)
unsigned int numFunctionCalls() const
void polint(std::vector< double >::iterator xArray, std::vector< double >::iterator yArray, double x, double &y, double &deltay) const
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const