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

StepDoublingRKStepper.cc
Go to the documentation of this file.
2 #include <stdexcept>
3 #include <cmath>
4 namespace Genfun {
5 
6 
7  StepDoublingRKStepper::StepDoublingRKStepper(const ButcherTableau & xtableau):tableau(xtableau) {
8  }
9 
11  }
12 
14  const RKIntegrator::RKData::Data & s,
16  std::vector<double> & errors) const {
17  const unsigned int nvar = s.variable.size();
18  RKIntegrator::RKData::Data d1(nvar),d2(nvar);
19 
20  doStep(data,s,d);
21  double dt = (d.time - s.time);
22  d1.time = s.time + dt/2.0;
23  d2.time = d.time;
24 
25  doStep(data, s,d1);
26  doStep(data,d1,d2);
27 
28  // Error estimate:
29  errors.resize(nvar);
30  for (size_t v=0;v<nvar;v++) errors[v]=fabs(d2.variable[v]-d.variable[v]);
31 
32  // Final correction:
33  for (size_t v=0;v<nvar;v++) d.variable[v] = d2.variable[v] + ((d2.variable[v]-d.variable[v])/double(std::pow(2.,(int)(tableau.order())-1)));
34 
35  }
36 
38  const RKIntegrator::RKData::Data & s,
39  RKIntegrator::RKData::Data & d) const {
40  // First step:
41  double h = (d.time - s.time);
42 
43 
44  if (h<=0) throw std::runtime_error ("SimpleRKStepper: negative stepsize");
45  const unsigned int nvar = s.variable.size();
46  // Compute all of the k's..:
47  //
48  std::vector<std::vector<double> >k(tableau.nSteps());
49  for (unsigned int i=0;i<tableau.nSteps();i++) {
50  k[i].resize(nvar,0);
51  Argument arg(nvar);
52  for (unsigned int v=0;v<nvar;v++) arg[v]=s.variable[v];
53  for (unsigned int j=0;j<i;j++) {
54  for (unsigned int v=0;v<nvar;v++) arg[v] += h*tableau.A(i,j)*k[j][v];
55  }
56  for (unsigned int v=0;v<nvar;v++) k[i][v]=(*data->_diffEqn[v])(arg);
57  }
58  //
59  // Final result.
60  //
61  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] = 0;
62  for (unsigned int i=0;i<tableau.nSteps();i++) {
63  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] += tableau.b(i)*k[i][v];
64  }
65  for (unsigned int v=0;v<nvar;v++) d.variable[v] =s.variable[v]+h*d.firstDerivative[v];
66 
67  }
68 
69 
70 
71 
73  return new StepDoublingRKStepper(*this);
74  }
75 
76  unsigned int StepDoublingRKStepper::order() const {
77  return tableau.order();
78  }
79 
80 
81 }