BALL
1.4.1
|
00001 // -*- Mode: C++; tab-width: 2; -*- 00002 // vi: set ts=2: 00003 // 00004 // $Id: rombergIntegrator.h,v 1.12 2003/08/26 08:04:22 oliver Exp $ 00005 // 00006 00007 #ifndef BALL_MATHS_ROMBERGINTEGRATOR_H 00008 #define BALL_MATHS_ROMBERGINTEGRATOR_H 00009 00010 #ifndef BALL_MATHS_NUMERICALINTERGRATOR_H 00011 # include <BALL/MATHS/numericalIntegrator.h> 00012 #endif 00013 00014 namespace BALL 00015 { 00019 template <typename Function, typename DataType> 00020 class RombergIntegrator: public NumericalIntegrator<Function, DataType> 00021 { 00022 public: 00023 00024 BALL_CREATE(RombergIntegrator) 00025 00026 00027 00028 00030 RombergIntegrator(float epsilon=1E-5, Size nsteps=1000); 00031 00033 RombergIntegrator(const RombergIntegrator& romint); 00034 00036 ~RombergIntegrator(); 00037 00039 00040 00041 00043 const RombergIntegrator& operator = (const RombergIntegrator& romint); 00044 00046 virtual void clear(); 00047 00049 void setEpsilon(float eps); 00050 00052 void setMaxNumSteps(Size mns); 00053 00055 00056 00057 00059 bool operator == (const RombergIntegrator& romint) const; 00060 00062 00063 00064 00070 DataType integrate(DataType from, DataType to); 00071 00072 00082 DataType trapezoid(DataType h, DataType from, DataType to); 00083 00085 00086 protected: 00087 00088 float epsilon_; 00089 Size maxNumSteps_; 00090 vector<DataType> result_; 00091 }; 00092 00093 template<typename Function, typename DataType> 00094 BALL_INLINE 00095 RombergIntegrator<Function, DataType>::RombergIntegrator(float eps, Size nsteps): NumericalIntegrator<Function, DataType>(), epsilon_(eps), maxNumSteps_(nsteps) 00096 { 00097 result_.reserve(maxNumSteps_ / 10); 00098 } 00099 00100 template<typename Function, typename DataType> 00101 BALL_INLINE 00102 RombergIntegrator<Function, DataType>::RombergIntegrator(const RombergIntegrator<Function, DataType>& romint):NumericalIntegrator<Function, DataType>(romint) 00103 { 00104 } 00105 00106 template<typename Function, typename DataType> 00107 BALL_INLINE 00108 RombergIntegrator<Function, DataType>::~RombergIntegrator() 00109 { 00110 } 00111 00112 template<typename Function, typename DataType> 00113 BALL_INLINE 00114 const RombergIntegrator<Function, DataType>& 00115 RombergIntegrator<Function, DataType>::operator = 00116 (const RombergIntegrator<Function, DataType>& romint) 00117 { 00118 function_ = romint.function_; 00119 maxNumSteps_ = romint.maxNumSteps_; 00120 epsilon_ = romint.epsilon_; 00121 result_ = romint.result_; 00122 return *this; 00123 } 00124 00125 template<typename Function, typename DataType> 00126 BALL_INLINE 00127 void RombergIntegrator<Function, DataType>::clear() 00128 { 00129 // Problem: function_.clear() might not exist... function_.clear(); 00130 } 00131 00132 template<typename Function, typename DataType> 00133 BALL_INLINE 00134 void RombergIntegrator<Function, DataType>::setEpsilon(float eps) 00135 { 00136 epsilon_ = eps; 00137 } 00138 00139 template<typename Function, typename DataType> 00140 BALL_INLINE 00141 void RombergIntegrator<Function, DataType>::setMaxNumSteps(Size nsteps) 00142 { 00143 maxNumSteps_ = nsteps; 00144 result_.reserve(maxNumSteps_ / 10); // we hope that we do not need more than 1/10 - th of the allowed operations 00145 } 00146 00147 template<typename Function, typename DataType> 00148 BALL_INLINE 00149 bool RombergIntegrator<Function, DataType>::operator == 00150 (const RombergIntegrator<Function, DataType> &romint) const 00151 00152 { 00153 return ((function_ == romint.function_) 00154 && (epsilon_ == romint.epsilon_ ) 00155 && (maxNumSteps_ == romint.maxNumSteps_) 00156 && (result_ == romint.result_ )); 00157 } 00158 00159 template<typename Function, typename DataType> 00160 BALL_INLINE 00161 DataType RombergIntegrator<Function, DataType>::trapezoid(DataType h, DataType from, DataType to) 00162 { 00163 DataType sum=0; 00164 DataType helper = (to - from); 00165 int count; 00166 00167 Size nsteps = (Size) (sqrt((helper*helper)/(h*h))); 00168 for (count=1; count<nsteps-1; count++) 00169 { 00170 sum +=function_(from+(count*h)); 00171 } 00172 00173 sum+=function_(from)+function_(to); 00174 sum*=h; 00175 00176 return sum; 00177 } 00178 00179 00180 template<typename Function, typename DataType> 00181 BALL_INLINE 00182 DataType RombergIntegrator<Function, DataType>::integrate 00183 (DataType from, DataType to) 00184 { 00185 float h = 1.; 00186 result_.push_back(trapezoid(h, from, to)); // this is the zeroth approximation 00187 00188 int i=1; 00189 int j=0; 00190 int count = 0; 00191 DataType dev; 00192 00193 do 00194 { 00195 result_.push_back(trapezoid(h/((float) i+1), from, to)); 00196 00197 for (j=1; j <= i; j++) 00198 { 00199 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])); 00200 count++; 00201 }; 00202 i++; 00203 dev = result_[((i-2)*(i-1))/2+(i-2)] - result_[((i-1)*(i))/2+(i-1)]; 00204 } while ( (sqrt(dev*dev) > epsilon_) && (count < maxNumSteps_)); 00205 00206 return (result_[((i-1)*(i))/2 + (i-1)]); 00207 } 00208 } // namespace BALL 00209 00210 #endif // BALL_MATHS_ROMBERGINTEGRATOR_H