CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

GenericFunctions/RKIntegrator.hh
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
3 //---------------------Runge-Kutte Integrator-------------------------------//
4 // //
5 // Class RKIntegrator //
6 // Joe Boudreau, November 2002 //
7 // //
8 // This is a Runge-Kutta Numerical integrator for a set of N autonomous //
9 // first order differential equations in N variables. The point is to //
10 // create one or more functions which are defined by A) the differential //
11 // equations governing their time evolution, and B) their values at time //
12 // t=0. //
13 // //
14 // You add differential eqns one at a time to this integrator. Each one //
15 // is a GENFUNCTION governing the time evolution of the i^th variable, and //
16 // should depend on all of the N variables, but not on the time //
17 // explicitly. You should add N differential equations in all. Each //
18 // time you add a differential equation the integrator creates a parameter //
19 // for you representing the starting value of the variable, and returns a //
20 // pointer. You may either set the values of that parameter to desired //
21 // values or else connect it to an external parameter if you wish to vary //
22 // the shape of the function by adjusting starting values. //
23 // //
24 // In addition, you may request the integrator to create a control //
25 // parameter. The control parameter may also be set, or connected. //
26 // It can be used in the equations that define the time evolution of the //
27 // variables. //
28 //--------------------------------------------------------------------------//
29 #ifndef RKIntegrator_h
30 #define RKIntegrator_h 1
34 #include <vector>
35 #include <set>
36 namespace Genfun {
37 
43  class RKIntegrator {
44 
45  public:
46 
47  // Some helper classes:
48  class RKFunction;
49  class RKData;
50  class RKStepper;
51 
52  // Constructor
53  RKIntegrator(const RKStepper *stepper=NULL);
54 
55  // Destructor
56  virtual ~RKIntegrator();
57 
58  // Add a differential equation governing the time evolution of the next variable.
59  // Get back a parameter representing the starting value of that variable. You
60  // can either arrange for that parameter to have the right starting value, or you
61  // can connect it to another parameter so that you may change it.
62  Parameter * addDiffEquation (const AbsFunction * diffEquation,
63  const std::string & variableName="anon",
64  double defStartingValue=0.0,
65  double startingValueMin=0.0,
66  double startingValueMax=0.0);
67 
68 
69  // Create a control parameter. You can then connnect this to some other
70  // parameter.
71  Parameter *createControlParameter (const std::string & variableName="anon",
72  double defStartingValue=0.0,
73  double startingValueMin=0.0,
74  double startingValueMax=0.0);
75 
76  // Get back a function. This function will now actually change as parameters
77  // are changed; this includes both control parameters and starting value
78  // parameters.
79  const RKFunction *getFunction(unsigned int i) const;
80 
81 
82  private:
83 
84  // It is illegal to assign an RKIntegrator
85  const RKIntegrator & operator=(const RKIntegrator &right);
86 
87  // It is illegal to copy an RKIntegrator
88  RKIntegrator(const RKIntegrator &right);
89 
90  // Here is the data, it belongs to the integrator and to the
91  // functions, and is reference counted:
92  RKData *_data;
93 
94 
95  // Here are the functions:
96  std::vector<const RKFunction *> _fcn;
97 
98 
99  };
100 
101 
102  class RKIntegrator::RKData : public Genfun::RCBase {
103 
104 
105  public:
106 
107  // Information about solution at each mesh point.
108  struct Data{
109 
110  std::vector<double> variable; // Solution
111  mutable std::vector<double> firstDerivative; // It's first derivative
112  double time; // time
113 
114  Data(int size): variable(size), firstDerivative(size), time(0) {}
115  bool operator < (const Data & right) const { return time < right.time; }
116  bool operator == (const Data & right) const { return time==right.time; }
117  };
118 
119  RKData();
120  void lock();
121  void recache();
122 
123  std::vector<Parameter *> _startingValParameter;
124  std::vector<double> _startingValParameterCache;
125 
126  std::vector <Parameter *> _controlParameter;
127  std::vector <double> _controlParameterCache;
128 
129  std::vector<const AbsFunction *> _diffEqn;
130  std::set<Data > _fx;
131  bool _locked;
132  const RKStepper *_stepper;
133  private:
134 
135  ~RKData();
136  friend class ImaginaryFriend; // Silence compiler warnings.
137 
138  };
139 
140  class RKIntegrator::RKFunction : public AbsFunction {
141 
143 
144  public:
145 
146  // Constructor
147  RKFunction(RKData *data, unsigned int index);
148 
149  // Destructor
150  virtual ~RKFunction();
151 
152  // Copy constructor
153  RKFunction(const RKFunction &right);
154 
155  // Retreive function value
156  virtual double operator ()(double argument) const;
157  virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
158 
159  private:
160 
161  // It is illegal to assign a RKFunction
162  const RKFunction & operator=(const RKFunction &right);
163 
164  // The shared data:
165  RKData *_data;
166  const unsigned int _index;
167 
168 };
169 
170 
171  // An abstract base class for steppers:
173  public:
174 
175  virtual ~RKStepper();
176  virtual void step (const RKIntegrator::RKData *data,
177  const RKIntegrator::RKData::Data & sdata,
179  double timeLimit=0) const =0;
180  virtual RKStepper *clone() const=0;
181 
182  };
183 
184 } // namespace Genfun
185 
186 #endif