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>
4namespace Genfun {
5 SimpleRKStepper::SimpleRKStepper(const ButcherTableau & mtableau,double xstepsize):
6 tableau(mtableau),
7 stepsize(xstepsize)
8 {
9 }
10
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}
unsigned int nSteps() const
double & A(unsigned int i, unsigned int j)
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