BALL 1.5.0
rombergIntegrator.h
Go to the documentation of this file.
1// -*- Mode: C++; tab-width: 2; -*-
2// vi: set ts=2:
3//
4// $Id: rombergIntegrator.h,v 1.12 2003/08/26 08:04:22 oliver Exp $
5//
6
7#ifndef BALL_MATHS_ROMBERGINTEGRATOR_H
8#define BALL_MATHS_ROMBERGINTEGRATOR_H
9
10#ifndef BALL_MATHS_NUMERICALINTERGRATOR_H
12#endif
13
14namespace BALL
15{
19 template <typename Function, typename DataType>
20 class RombergIntegrator: public NumericalIntegrator<Function, DataType>
21 {
22 public:
23
25
26
28
30 RombergIntegrator(float epsilon=1E-5, Size nsteps=1000);
31
34
37
39
41
44
46 virtual void clear();
47
49 void setEpsilon(float eps);
50
52 void setMaxNumSteps(Size mns);
53
55
57
59 bool operator == (const RombergIntegrator& romint) const;
60
62
64
70 DataType integrate(DataType from, DataType to);
71
72
82 DataType trapezoid(DataType h, DataType from, DataType to);
83
85
86 protected:
87
88 float epsilon_;
90 vector<DataType> result_;
91 };
92
93 template<typename Function, typename DataType>
95 RombergIntegrator<Function, DataType>::RombergIntegrator(float eps, Size nsteps): NumericalIntegrator<Function, DataType>(), epsilon_(eps), maxNumSteps_(nsteps)
96 {
97 result_.reserve(maxNumSteps_ / 10);
98 }
99
100 template<typename Function, typename DataType>
103 {
104 }
105
106 template<typename Function, typename DataType>
109 {
110 }
111
112 template<typename Function, typename DataType>
117 {
118 function_ = romint.function_;
119 maxNumSteps_ = romint.maxNumSteps_;
120 epsilon_ = romint.epsilon_;
121 result_ = romint.result_;
122 return *this;
123 }
124
125 template<typename Function, typename DataType>
128 {
129 // Problem: function_.clear() might not exist... function_.clear();
130 }
131
132 template<typename Function, typename DataType>
135 {
136 epsilon_ = eps;
137 }
138
139 template<typename Function, typename DataType>
142 {
143 maxNumSteps_ = nsteps;
144 result_.reserve(maxNumSteps_ / 10); // we hope that we do not need more than 1/10 - th of the allowed operations
145 }
146
147 template<typename Function, typename DataType>
150 (const RombergIntegrator<Function, DataType> &romint) const
151
152 {
153 return ((function_ == romint.function_)
154 && (epsilon_ == romint.epsilon_ )
155 && (maxNumSteps_ == romint.maxNumSteps_)
156 && (result_ == romint.result_ ));
157 }
158
159 template<typename Function, typename DataType>
161 DataType RombergIntegrator<Function, DataType>::trapezoid(DataType h, DataType from, DataType to)
162 {
163 DataType sum=0;
164 DataType helper = (to - from);
165 int count;
166
167 Size nsteps = (Size) (sqrt((helper*helper)/(h*h)));
168 for (count=1; count<nsteps-1; count++)
169 {
170 sum +=function_(from+(count*h));
171 }
172
173 sum+=function_(from)+function_(to);
174 sum*=h;
175
176 return sum;
177 }
178
179
180 template<typename Function, typename DataType>
183 (DataType from, DataType to)
184 {
185 float h = 1.;
186 result_.push_back(trapezoid(h, from, to)); // this is the zeroth approximation
187
188 int i=1;
189 int j=0;
190 int count = 0;
191 DataType dev;
192
193 do
194 {
195 result_.push_back(trapezoid(h/((float) i+1), from, to));
196
197 for (j=1; j <= i; j++)
198 {
199 result_.push_back(result_[(i*(i+1))/2 + (j-1)] + 1. / (pow(4, j) - 1) * (result_[(i*(i+1))/2 + j-1 - result_[((i-1)*i)/2+j-1]));
200 count++;
201 };
202 i++;
203 dev = result_[((i-2)*(i-1))/2+(i-2)] - result_[((i-1)*(i))/2+(i-1)];
204 } while ( (sqrt(dev*dev) > epsilon_) && (count < maxNumSteps_));
205
206 return (result_[((i-1)*(i))/2 + (i-1)]);
207 }
208} // namespace BALL
209
210#endif // BALL_MATHS_ROMBERGINTEGRATOR_H
#define BALL_CREATE(name)
Definition: create.h:62
#define BALL_INLINE
Definition: debug.h:15
Definition: constants.h:13
BALL_SIZE_TYPE Size
BALL_EXTERN_VARIABLE const double h
Definition: constants.h:102
BALL_EXTERN_VARIABLE const double E
Euler's number - base of the natural logarithm.
Definition: constants.h:38
void setMaxNumSteps(Size mns)
Set the maximum number of steps we want to use in computation.
DataType integrate(DataType from, DataType to)
RombergIntegrator(float epsilon=1E-5, Size nsteps=1000)
Default constructor.
const RombergIntegrator & operator=(const RombergIntegrator &romint)
Assignment operator.
DataType trapezoid(DataType h, DataType from, DataType to)
bool operator==(const RombergIntegrator &romint) const
Equality operator.
virtual void clear()
Clear method.
vector< DataType > result_
void setEpsilon(float eps)
Set the upper bound for the error we want to allow.