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

EmbeddedRKStepper.cc
Go to the documentation of this file.
3 #include <stdexcept>
4 namespace Genfun {
5 
6 
8  tableau(mtableau){
9  }
10 
12  }
13 
15  const RKIntegrator::RKData::Data & s,
17  std::vector<double> & errors) const {
18 
19  // First step:
20  double h = d.time - s.time;
21  if (h<=0) throw std::runtime_error ("Runtime error in RKIntegrator (zero or negative stepsize)");
22  unsigned int nvar = s.variable.size();
23 
24  // Compute all of the k's..:
25  //
26  std::vector<std::vector<double> >k(tableau.nSteps());
27  for (unsigned int i=0;i<tableau.nSteps();i++) {
28  k[i].resize(nvar,0);
29  Argument arg(nvar);
30  for (unsigned int v=0;v<nvar;v++) arg[v]=s.variable[v];
31  for (unsigned int j=0;j<i;j++) {
32  for (unsigned int v=0;v<nvar;v++) arg[v] += h*tableau.A(i,j)*k[j][v];
33  }
34  for (unsigned int v=0;v<nvar;v++) k[i][v]=(*data->_diffEqn[v])(arg);
35  }
36  //
37  // Final result.
38  //
39  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] = 0;
40  for (unsigned int i=0;i<tableau.nSteps();i++) {
41  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] += tableau.b(i)*k[i][v];
42  }
43  for (unsigned int v=0;v<nvar;v++) d.variable[v] =s.variable[v]+h*d.firstDerivative[v];
44  //
45  //
46  //
47  errors.resize(nvar);
48  for (unsigned int v=0;v<nvar;v++) errors[v] = 0;
49  for (unsigned int i=0;i<tableau.nSteps();i++) {
50  for (unsigned int v=0;v<nvar;v++) errors[v] += (h*(tableau.bHat(i)-tableau.b(i))*k[i][v]);
51  }
52  return;
53  }
54 
56  return new EmbeddedRKStepper(*this);
57  }
58 
59  unsigned int EmbeddedRKStepper::order() const {
60  return tableau.order();
61  }
62 }
double & A(unsigned int i, unsigned int j)
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, std::vector< double > &errors) const
virtual EmbeddedRKStepper * clone() const
double & bHat(unsigned int i)
virtual unsigned int order() const
unsigned int order() const
unsigned int nSteps() const
std::vector< const AbsFunction * > _diffEqn
EmbeddedRKStepper(const ExtendedButcherTableau &tableau=CashKarpXtTableau())
double & b(unsigned int i)