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

RKIntegrator.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
5 #include <climits>
6 #include <cmath> // for pow()
7 #include <stdexcept>
8 namespace Genfun {
9 FUNCTION_OBJECT_IMP(RKIntegrator::RKFunction)
10 
11 RKIntegrator::RKFunction::RKFunction(RKData *data, unsigned int index)
12  :_data(data),
13  _index(index)
14 {
15  _data->ref();
16 }
17 
19 {
20  _data->unref();
21 }
22 
24  :_data(right._data),
25  _index(right._index)
26 {
27  _data->ref();
28 }
29 
30 
31 double RKIntegrator::RKFunction::operator() (double t) const {
32  if (t<0) return 0;
33  if (!_data->_locked) _data->lock();
34 
35 
36  // Do this first, thereafter, just read the cache
37  _data->recache();
38 
39 
40  // If the cache is empty, make an entry for t=0;
41  size_t nvar = _data->_startingValParameter.size();
42  if (_data->_fx.empty()) {
43  RKData::Data d(nvar);
44  d.time=0;
45  Argument arg(nvar);
46  for (size_t f=0;f<nvar;f++) {
48  arg[f]=d.variable[f];
49  }
50  _data->_fx.insert(d);
51  }
52 
53 
54  RKData::Data dt(nvar);
55  dt.time=t;
56  std::set<RKData::Data>::iterator s=_data->_fx.find(dt);
57  if (s!=_data->_fx.end()) {
58  // Then, there is nothing to do. Don't touch the
59  // list. Just get the variable:
60  return (*s).variable[_index];
61  }
62  else {
63  s=_data->_fx.lower_bound(dt);
64  // Back up:
65  if (!(s!=_data->_fx.begin())) throw std::runtime_error("Runtime error in RKIntegrator");
66  s--;
67 
68  //std::vector<double> errors;
69  rkstep(*s, dt);
70  _data->_fx.insert(s,dt);
71 
72  return dt.variable[_index];
73  }
74 }
75 
76 
77 RKIntegrator::RKData::RKData():_locked(false) {
78 }
79 
80 RKIntegrator::RKData::~RKData() {
81  for (size_t i=0;i<_startingValParameter.size();i++) delete _startingValParameter[i];
82  for (size_t i=0;i<_controlParameter.size();i++) delete _controlParameter[i];
83  for (size_t i=0;i<_diffEqn.size(); i++) delete _diffEqn[i];
84 }
85 
87  _data(new RKData())
88 {
89  _data->ref();
90 }
91 
93  _data->unref();
94  for (size_t i=0;i<_fcn.size();i++) delete _fcn[i];
95 }
96 
98  const std::string &variableName,
99  double defStartingValue,
100  double defValueMin,
101  double defValueMax) {
102  Parameter *par = new Parameter(variableName, defStartingValue, defValueMin, defValueMax);
103  _data->_startingValParameter.push_back(par);
104  _data->_diffEqn.push_back(diffEquation->clone());
105  _data->_startingValParameterCache.push_back(defStartingValue);
106  _fcn.push_back(new RKFunction(_data,_fcn.size()));
107  return par;
108 }
109 
110 
111 
112 
113 
114 Parameter * RKIntegrator::createControlParameter (const std::string & variableName,
115  double defStartingValue,
116  double startingValueMin,
117  double startingValueMax) {
118 
119  Parameter *par = new Parameter(variableName, defStartingValue, startingValueMin, startingValueMax);
120  _data->_controlParameter.push_back(par);
121  _data->_controlParameterCache.push_back(defStartingValue);
122  return par;
123 
124 }
125 
126 
127 
129  return _fcn[i];
130 }
131 
132 
133 
135  if (!_locked) {
136  unsigned int size = _diffEqn.size();
137  for (size_t i=0;i<size;i++) {
138  if (!(_diffEqn[i]->dimensionality()==size)) throw std::runtime_error("Runtime error in RKIntegrator");
139  }
140  _locked=true;
141  }
142 }
143 
144 
145 
147 
148  bool stale=false;
149  if (!stale) {
150  for (size_t p=0;p<_startingValParameter.size();p++) {
151  if (_startingValParameter[p]->getValue()!=_startingValParameterCache[p]) {
152  _startingValParameterCache[p]=_startingValParameter[p]->getValue();
153  stale=true;
154  break;
155  }
156  }
157  }
158 
159  if (!stale) {
160  for (size_t p=0;p<_controlParameter.size();p++) {
161  if (_controlParameter[p]->getValue()!=_controlParameterCache[p]) {
162  _controlParameterCache[p]=_controlParameter[p]->getValue();
163  stale=true;
164  break;
165  }
166  }
167  }
168 
169  if (stale) {
170  _fx.erase(_fx.begin(),_fx.end());
171  }
172 
173 }
174 
175 
176 
177 void RKIntegrator::RKFunction::rkstep(const RKIntegrator::RKData::Data & s, RKIntegrator::RKData::Data & d) const {
178  //
179  // Adaptive stepsize control
180  //
181  const int nvar = s.variable.size();
182  const double eps = 1.0E-6;
183  const double SAFETY = 0.9;
184  const double PSHRNK = -0.25;
185  const double PGROW = -0.20;
186  const double ERRCON = -1.89E-4;
187  const double TINY = 1.0E-30;
188  double hnext;
189 
190 
191  RKData::Data Tmp0(nvar),Tmp1(nvar);
192  Tmp0=s;
193  Tmp1=d;
194  bool done=false;
195  while (1) { // "Same time as"...
196 
197  //--------------------------------------//
198  // Take one step, from Tmp0 to Tmp1: //
199  //--------------------------------------//
200 
201  double h = Tmp1.time - Tmp0.time;
202  while (1) {
203  std::vector<double> errors;
204 
205  rkck(Tmp0, Tmp1, errors);
206  for (size_t e=0;e<errors.size();e++) {
207  errors[e] = fabs(errors[e]) / (fabs(Tmp0.variable[e]) + fabs(Tmp0.firstDerivative[e]*h) + TINY);
208  }
209  double emax = (*std::max_element(errors.begin(),errors.end()))/eps;
210  if (emax > 1) {
211  h = std::max(SAFETY*h*pow(emax,PSHRNK),0.1*h);
212  if (!(((float) Tmp0.time+h - (float) Tmp0.time) > 0) ) {
213  std::cerr << "Warning, RK Integrator step underflow" << std::endl;
214  }
215  Tmp1.time = Tmp0.time+h;
216  continue;
217  }
218  else {
219 
220  if (emax > ERRCON) {
221  hnext = SAFETY*h*pow(emax,PGROW);
222  }
223  else {
224  hnext = 5.0*h;
225  }
226  if (Tmp1==d) {
227  done=true;
228  break;
229  }
230  else {
231  Tmp0=Tmp1;
232  Tmp1.time = std::min(Tmp0.time + hnext, d.time);
233  break;
234  }
235  }
236  }
237 
238  //--------------------------------------//
239  // End of Step. //
240  //--------------------------------------//
241 
242  if (done) break;
243  }
244  d=Tmp1;
245 }
246 
247 void RKIntegrator::RKFunction::rkck(const RKIntegrator::RKData::Data & s, RKIntegrator::RKData::Data & d, std::vector<double> & errors) const {
248 
249 #ifdef NONAUTONOMOUS_EQUATIONS
250  static const double
251  a2=0.2,
252  a3=0.3,
253  a4=0.6,
254  a5=1.0,
255  a6=0.875;
256 #endif
257 
258  static const double
259  b21=0.2,
260  b31=3.0/40.0,
261  b32=9.0/40.0,
262  b41=0.3,
263  b42=-0.9,
264  b43=1.2,
265  b51=-11.0/54.0,
266  b52=2.5,
267  b53=-70.0/27.0,
268  b54=35.0/27.0,
269  b61=1631.0/55296.0,
270  b62=175.0/512.0,
271  b63=575.0/13824.0,
272  b64=44275.0/110592.0,
273  b65=253.0/4096.0;
274 
275  static const double
276  c1=37.0/378.0,
277  c3=250.0/621.0,
278  c4=125.0/594.0,
279  c6=512.0/1771.0;
280 
281  static const double
282  dc1=c1-2825.0/27648.0,
283  dc3=c3-18575.0/48384.0,
284  dc4=c4-13525.0/55296.0,
285  dc5=-277.0/14336.0,
286  dc6=c6 - 0.25;
287 
288  // First step:
289  double h = d.time - s.time;
290  if (h<=0) throw std::runtime_error ("Runtime error in RKIntegrator (zero or negative stepsize)");
291  unsigned int nv = s.variable.size();
292  Argument arg(nv), arg0(nv), d1(nv),d2(nv), d3(nv), d4(nv), d5(nv), d6(nv);
293 
294 
295  for (size_t v=0;v<nv;v++) { arg0[v]=s.variable[v];}
296 
297  if (!s.dcalc) {
298  for (size_t v=0;v<nv;v++) {d1[v]=(*_data->_diffEqn[v])(arg0);}
299  for (size_t v=0;v<nv;v++) {s.firstDerivative[v]=d1[v];}
300  s.dcalc=true;
301  }
302  else {
303  for (size_t v=0;v<nv;v++) { d1[v] = s.firstDerivative[v];}
304  }
305 
306 
307  for (size_t v=0;v<nv;v++) { arg[v] = arg0[v] + b21*h*d1[v];}
308 
309  for (size_t v=0;v<nv;v++) { d2[v] = (*_data->_diffEqn[v])(arg);}
310  for (size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b31*d1[v]+b32*d2[v]);}
311 
312 
313  for (size_t v=0;v<nv;v++) { d3[v] = (*_data->_diffEqn[v])(arg);}
314  for (size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b41*d1[v]+b42*d2[v]+b43*d3[v]);}
315 
316  for (size_t v=0;v<nv;v++) { d4[v] = (*_data->_diffEqn[v])(arg);}
317  for (size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b51*d1[v]+b52*d2[v]+b53*d3[v] + b54*d4[v]);}
318 
319  for (size_t v=0;v<nv;v++) { d5[v] = (*_data->_diffEqn[v])(arg);}
320  for (size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b61*d1[v]+b62*d2[v]+b63*d3[v] + b64*d4[v] + b65*d5[v]);}
321 
322  for (size_t v=0;v<nv;v++) { d6[v] = (*_data->_diffEqn[v])(arg);}
323 
324  for (size_t v=0;v<nv;v++) { d.variable[v] = arg0[v] + h*(c1*d1[v]+c3*d3[v]+c4*d4[v]+c6*d6[v]);}
325  errors.erase(errors.begin(),errors.end());
326 
327  for (size_t v=0;v<nv;v++) { errors.push_back(h*(dc1*d1[v]+dc3*d3[v]+dc4*d4[v]+dc5*d5[v]+dc6*d6[v]));}
328 }
329 
330 
331 
332 } // namespace Genfun