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
11namespace 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
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);
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
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
virtual unsigned int stepMultiplier() const =0
virtual unsigned int numFunctionCalls() const =0
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const =0
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const
virtual double integrate(const AbsFunction &function, double a, double b, unsigned int j) const
void polint(std::vector< double >::iterator xArray, std::vector< double >::iterator yArray, double x, double &y, double &deltay) const
DefiniteIntegral(double a, double b, Type=CLOSED)
void setMinOrder(unsigned int order)
virtual double operator[](const AbsFunction &function) const
unsigned int numFunctionCalls() const
DefiniteIntegral & operator=(const DefiniteIntegral &)
void setMaxIter(unsigned int maxIter)