CLHEP VERSION 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 class RKStepper; 00051 00052 // Constructor 00053 RKIntegrator(const RKStepper *stepper=NULL); 00054 00055 // Destructor 00056 virtual ~RKIntegrator(); 00057 00058 // Add a differential equation governing the time evolution of the next variable. 00059 // Get back a parameter representing the starting value of that variable. You 00060 // can either arrange for that parameter to have the right starting value, or you 00061 // can connect it to another parameter so that you may change it. 00062 Parameter * addDiffEquation (const AbsFunction * diffEquation, 00063 const std::string & variableName="anon", 00064 double defStartingValue=0.0, 00065 double startingValueMin=0.0, 00066 double startingValueMax=0.0); 00067 00068 00069 // Create a control parameter. You can then connnect this to some other 00070 // parameter. 00071 Parameter *createControlParameter (const std::string & variableName="anon", 00072 double defStartingValue=0.0, 00073 double startingValueMin=0.0, 00074 double startingValueMax=0.0); 00075 00076 // Get back a function. This function will now actually change as parameters 00077 // are changed; this includes both control parameters and starting value 00078 // parameters. 00079 const RKFunction *getFunction(unsigned int i) const; 00080 00081 00082 private: 00083 00084 // It is illegal to assign an RKIntegrator 00085 const RKIntegrator & operator=(const RKIntegrator &right); 00086 00087 // It is illegal to copy an RKIntegrator 00088 RKIntegrator(const RKIntegrator &right); 00089 00090 // Here is the data, it belongs to the integrator and to the 00091 // functions, and is reference counted: 00092 RKData *_data; 00093 00094 00095 // Here are the functions: 00096 std::vector<const RKFunction *> _fcn; 00097 00098 00099 }; 00100 00101 00102 class RKIntegrator::RKData : public Genfun::RCBase { 00103 00104 00105 public: 00106 00107 // Information about solution at each mesh point. 00108 struct Data{ 00109 00110 std::vector<double> variable; // Solution 00111 mutable std::vector<double> firstDerivative; // It's first derivative 00112 double time; // time 00113 00114 Data(int size): variable(size), firstDerivative(size), time(0) {} 00115 bool operator < (const Data & right) const { return time < right.time; } 00116 bool operator == (const Data & right) const { return time==right.time; } 00117 }; 00118 00119 RKData(); 00120 void lock(); 00121 void recache(); 00122 00123 std::vector<Parameter *> _startingValParameter; 00124 std::vector<double> _startingValParameterCache; 00125 00126 std::vector <Parameter *> _controlParameter; 00127 std::vector <double> _controlParameterCache; 00128 00129 std::vector<const AbsFunction *> _diffEqn; 00130 std::set<Data > _fx; 00131 bool _locked; 00132 const RKStepper *_stepper; 00133 private: 00134 00135 ~RKData(); 00136 friend class ImaginaryFriend; // Silence compiler warnings. 00137 00138 }; 00139 00140 class RKIntegrator::RKFunction : public AbsFunction { 00141 00142 FUNCTION_OBJECT_DEF(RKFunction) 00143 00144 public: 00145 00146 // Constructor 00147 RKFunction(RKData *data, unsigned int index); 00148 00149 // Destructor 00150 virtual ~RKFunction(); 00151 00152 // Copy constructor 00153 RKFunction(const RKFunction &right); 00154 00155 // Retreive function value 00156 virtual double operator ()(double argument) const; 00157 virtual double operator ()(const Argument & a) const {return operator() (a[0]);} 00158 00159 private: 00160 00161 // It is illegal to assign a RKFunction 00162 const RKFunction & operator=(const RKFunction &right); 00163 00164 // The shared data: 00165 RKData *_data; 00166 const unsigned int _index; 00167 00168 }; 00169 00170 00171 // An abstract base class for steppers: 00172 class RKIntegrator::RKStepper { 00173 public: 00174 00175 virtual ~RKStepper(); 00176 virtual void step (const RKIntegrator::RKData *data, 00177 const RKIntegrator::RKData::Data & sdata, 00178 RKIntegrator::RKData::Data & ddata, 00179 double timeLimit=0) const =0; 00180 virtual RKStepper *clone() const=0; 00181 00182 }; 00183 00184 } // namespace Genfun 00185 00186 #endif