CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: 00003 //---------------------Runge-Kutte Integrator-------------------------------// 00004 // // 00005 // Class RKIntegrator // 00006 // Joe Boudreau, November 2002 // 00007 // // 00008 // This is a Runge-Kutta Numerical integrator for a set of N autonomous // 00009 // first order differential equations in N variables. The point is to // 00010 // create one or more functions which are defined by A) the differential // 00011 // equations governing their time evolution, and B) their values at time // 00012 // t=0. // 00013 // // 00014 // You add differential eqns one at a time to this integrator. Each one // 00015 // is a GENFUNCTION governing the time evolution of the i^th variable, and // 00016 // should depend on all of the N variables, but not on the time // 00017 // explicitly. You should add N differential equations in all. Each // 00018 // time you add a differential equation the integrator creates a parameter // 00019 // for you representing the starting value of the variable, and returns a // 00020 // pointer. You may either set the values of that parameter to desired // 00021 // values or else connect it to an external parameter if you wish to vary // 00022 // the shape of the function by adjusting starting values. // 00023 // // 00024 // In addition, you may request the integrator to create a control // 00025 // parameter. The control parameter may also be set, or connected. // 00026 // It can be used in the equations that define the time evolution of the // 00027 // variables. // 00028 //--------------------------------------------------------------------------// 00029 #ifndef RKIntegrator_h 00030 #define RKIntegrator_h 1 00031 #include "CLHEP/GenericFunctions/AbsFunction.hh" 00032 #include "CLHEP/GenericFunctions/Parameter.hh" 00033 #include "CLHEP/GenericFunctions/RCBase.hh" 00034 #include <vector> 00035 #include <set> 00036 namespace Genfun { 00037 00043 class RKIntegrator { 00044 00045 public: 00046 00047 // Some helper classes: 00048 class RKFunction; 00049 class RKData; 00050 00051 // Constructor 00052 RKIntegrator(); 00053 00054 // Destructor 00055 virtual ~RKIntegrator(); 00056 00057 // Add a differential equation governing the time evolution of the next variable. 00058 // Get back a parameter representing the starting value of that variable. You 00059 // can either arrange for that parameter to have the right starting value, or you 00060 // can connect it to another parameter so that you may change it. 00061 Parameter * addDiffEquation (const AbsFunction * diffEquation, 00062 const std::string & variableName="anon", 00063 double defStartingValue=0.0, 00064 double startingValueMin=0.0, 00065 double startingValueMax=0.0); 00066 00067 00068 // Create a control parameter. You can then connnect this to some other 00069 // parameter. 00070 Parameter *createControlParameter (const std::string & variableName="anon", 00071 double defStartingValue=0.0, 00072 double startingValueMin=0.0, 00073 double startingValueMax=0.0); 00074 00075 // Get back a function. This function will now actually change as parameters 00076 // are changed; this includes both control parameters and starting value 00077 // parameters. 00078 const RKFunction *getFunction(unsigned int i) const; 00079 00080 00081 private: 00082 00083 // It is illegal to assign an adjustable constant 00084 const RKIntegrator & operator=(const RKIntegrator &right); 00085 00086 // It is illegal to copy an RKIntegrator 00087 RKIntegrator(const RKIntegrator &right); 00088 00089 // Here is the data, it belongs to the integrator and to the 00090 // functions, and is reference counted: 00091 RKData *_data; 00092 00093 00094 // Here are the functions: 00095 std::vector<const RKFunction *> _fcn; 00096 00097 00098 }; 00099 00100 00101 class RKIntegrator::RKData : public Genfun::RCBase { 00102 00103 00104 public: 00105 00106 struct Data{ 00107 00108 std::vector<double> variable; 00109 mutable std::vector<double> firstDerivative; 00110 double time; 00111 mutable bool dcalc; 00112 00113 Data(int size): variable(size), firstDerivative(size), time(0), dcalc(false) {} 00114 bool operator < (const Data & right) const { return time < right.time; } 00115 bool operator == (const Data & right) const { return time==right.time; } 00116 }; 00117 00118 RKData(); 00119 void lock(); 00120 void recache(); 00121 00122 std::vector<Parameter *> _startingValParameter; 00123 std::vector<double> _startingValParameterCache; 00124 00125 std::vector <Parameter *> _controlParameter; 00126 std::vector <double> _controlParameterCache; 00127 00128 std::vector<const AbsFunction *> _diffEqn; 00129 std::set<Data > _fx; 00130 bool _locked; 00131 00132 private: 00133 00134 ~RKData(); 00135 friend class ImaginaryFriend; // Silence compiler warnings. 00136 00137 }; 00138 00139 class RKIntegrator::RKFunction : public AbsFunction { 00140 00141 FUNCTION_OBJECT_DEF(RKFunction) 00142 00143 public: 00144 00145 // Constructor 00146 RKFunction(RKData *data, unsigned int index); 00147 00148 // Destructor 00149 virtual ~RKFunction(); 00150 00151 // Copy constructor 00152 RKFunction(const RKFunction &right); 00153 00154 // Retreive function value 00155 virtual double operator ()(double argument) const; 00156 virtual double operator ()(const Argument & a) const {return operator() (a[0]);} 00157 00158 private: 00159 00160 // It is illegal to assign a RKFunction 00161 const RKFunction & operator=(const RKFunction &right); 00162 00163 // The shared data: 00164 RKData *_data; 00165 const unsigned int _index; 00166 00167 void rk4 (const RKData::Data & sdata, RKData::Data & ddata) const; 00168 00169 void rkstep (const RKData::Data & sdata, RKData::Data & ddata) const; 00170 void rkck (const RKData::Data & sdata, RKData::Data & ddata, std::vector<double> & errors) const; 00171 }; 00172 00173 00174 } // namespace Genfun 00175 00176 #endif