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 }
Genfun::RKIntegrator::RKData::Data::time
double time
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:112
Genfun::ButcherTableau::b
double & b(unsigned int i)
Genfun::StepDoublingRKStepper::order
virtual unsigned int order() const
Definition: StepDoublingRKStepper.cc:76
Genfun::ButcherTableau::A
double & A(unsigned int i, unsigned int j)
StepDoublingRKStepper.hh
Genfun::StepDoublingRKStepper::step
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, std::vector< double > &errors) const
Definition: StepDoublingRKStepper.cc:13
Genfun::StepDoublingRKStepper::~StepDoublingRKStepper
virtual ~StepDoublingRKStepper()
Definition: StepDoublingRKStepper.cc:10
Genfun::RKIntegrator::RKData::Data
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:108
Genfun::RKIntegrator::RKData::Data::variable
std::vector< double > variable
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:110
Genfun::StepDoublingRKStepper
Definition: CLHEP/GenericFunctions/StepDoublingRKStepper.hh:11
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
Genfun::RKIntegrator::RKData
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:102
Genfun::ButcherTableau::order
unsigned int order() const
Genfun::Argument
Definition: CLHEP/GenericFunctions/Argument.hh:17
j
long j
Definition: JamesRandomSeeding.txt:28
Genfun::RKIntegrator::RKData::_diffEqn
std::vector< const AbsFunction * > _diffEqn
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:129
Genfun::RKIntegrator::RKData::Data::firstDerivative
std::vector< double > firstDerivative
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:111
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
i
long i
Definition: JamesRandomSeeding.txt:27
Genfun::StepDoublingRKStepper::StepDoublingRKStepper
StepDoublingRKStepper(const ButcherTableau &tableau)
Definition: StepDoublingRKStepper.cc:7
Genfun::StepDoublingRKStepper::clone
virtual StepDoublingRKStepper * clone() const
Definition: StepDoublingRKStepper.cc:72
for
for(n=1 ;n< 98 ;n++)
Definition: JamesRandomSeeding.txt:34
Genfun::ButcherTableau::nSteps
unsigned int nSteps() const
k
long k
Definition: JamesRandomSeeding.txt:29
Genfun::ButcherTableau
Definition: CLHEP/GenericFunctions/ButcherTableau.hh:23
Genfun::StepDoublingRKStepper::doStep
void doStep(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &s, RKIntegrator::RKData::Data &d) const
Definition: StepDoublingRKStepper.cc:37
Genfun
Definition: CLHEP/GenericFunctions/Abs.hh:14