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

SimpleRKStepper.cc
Go to the documentation of this file.
2 #include <cmath>
3 #include <stdexcept>
4 namespace Genfun {
5  SimpleRKStepper::SimpleRKStepper(const ButcherTableau & mtableau,double xstepsize):
6  tableau(mtableau),
7  stepsize(xstepsize)
8  {
9  }
10 
12  const RKIntegrator::RKData::Data & s,
14  double timeLimit ) const {
15  const double h = timeLimit==0 ? stepsize : timeLimit - s.time;
16  if (h<=0) throw std::runtime_error ("SimpleRKStepper: negative stepsize");
17  const unsigned int nvar = s.variable.size();
18  // Compute all of the k's..:
19  //
20  std::vector<std::vector<double> >k(tableau.nSteps());
21  for (unsigned int i=0;i<tableau.nSteps();i++) {
22  k[i].resize(nvar,0);
23  Argument arg(nvar);
24  for (unsigned int v=0;v<nvar;v++) arg[v]=s.variable[v];
25  for (unsigned int j=0;j<i;j++) {
26  for (unsigned int v=0;v<nvar;v++) arg[v] += h*tableau.A(i,j)*k[j][v];
27  }
28  for (unsigned int v=0;v<nvar;v++) k[i][v]=(*data->_diffEqn[v])(arg);
29  }
30  //
31  // Final result.
32  //
33  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] = 0;
34  for (unsigned int i=0;i<tableau.nSteps();i++) {
35  for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] += tableau.b(i)*k[i][v];
36  }
37  for (unsigned int v=0;v<nvar;v++) d.variable[v] =s.variable[v]+h*d.firstDerivative[v];
38  d.time = timeLimit==0 ? s.time + h : timeLimit;
39 
40  }
41 
43 
45  return new SimpleRKStepper(*this);
46  }
47 }
double & A(unsigned int i, unsigned int j)
unsigned int nSteps() const
double & b(unsigned int i)
std::vector< const AbsFunction * > _diffEqn
SimpleRKStepper(const ButcherTableau &tableau, double stepsize)
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit) const
virtual SimpleRKStepper * clone() const