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 
14 namespace BALL
15 {
19  template <typename Function, typename DataType>
20  class RombergIntegrator: public NumericalIntegrator<Function, DataType>
21  {
22  public:
23 
25 
26 
27 
28 
30  RombergIntegrator(float epsilon=1E-5, Size nsteps=1000);
31 
33  RombergIntegrator(const RombergIntegrator& romint);
34 
37 
39 
41 
43  const RombergIntegrator& operator = (const RombergIntegrator& romint);
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
BALL::RombergIntegrator::trapezoid
DataType trapezoid(DataType h, DataType from, DataType to)
Definition: rombergIntegrator.h:161
BALL_INLINE
#define BALL_INLINE
Definition: debug.h:15
BALL::RombergIntegrator::setEpsilon
void setEpsilon(float eps)
Set the upper bound for the error we want to allow.
Definition: rombergIntegrator.h:134
BALL::RombergIntegrator
Definition: rombergIntegrator.h:21
BALL::RombergIntegrator::operator=
const RombergIntegrator & operator=(const RombergIntegrator &romint)
Assignment operator.
Definition: rombergIntegrator.h:116
BALL::RombergIntegrator::maxNumSteps_
Size maxNumSteps_
Definition: rombergIntegrator.h:89
BALL::Size
BALL_SIZE_TYPE Size
Definition: COMMON/global.h:114
BALL::RombergIntegrator::result_
vector< DataType > result_
Definition: rombergIntegrator.h:90
BALL::RombergIntegrator::operator==
bool operator==(const RombergIntegrator &romint) const
Equality operator.
Definition: rombergIntegrator.h:150
BALL
Definition: constants.h:13
numericalIntegrator.h
BALL::RombergIntegrator::integrate
DataType integrate(DataType from, DataType to)
Definition: rombergIntegrator.h:183
BALL::RombergIntegrator::setMaxNumSteps
void setMaxNumSteps(Size mns)
Set the maximum number of steps we want to use in computation.
Definition: rombergIntegrator.h:141
BALL::NumericalIntegrator::function_
Function function_
Definition: numericalIntegrator.h:97
BALL_SIZE_TYPE
BALL::RombergIntegrator::~RombergIntegrator
~RombergIntegrator()
Destructor.
Definition: rombergIntegrator.h:108
BALL::RombergIntegrator::RombergIntegrator
RombergIntegrator(float epsilon=1E-5, Size nsteps=1000)
Default constructor.
Definition: rombergIntegrator.h:95
BALL::NumericalIntegrator
Definition: numericalIntegrator.h:21
BALL::Constants::E
BALL_EXTERN_VARIABLE const double E
Euler's number - base of the natural logarithm.
Definition: constants.h:38
BALL_CREATE
#define BALL_CREATE(name)
Definition: create.h:62
BALL::Constants::h
BALL_EXTERN_VARIABLE const double h
Definition: constants.h:102
BALL::RombergIntegrator::clear
virtual void clear()
Clear method.
Definition: rombergIntegrator.h:127
BALL::RombergIntegrator::epsilon_
float epsilon_
Definition: rombergIntegrator.h:88