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

EmbeddedRKStepper.cc
Go to the documentation of this file.
3#include <stdexcept>
4namespace Genfun {
5
6
8 tableau(mtableau){
9 }
10
12 }
13
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}
EmbeddedRKStepper(const ExtendedButcherTableau &tableau=CashKarpXtTableau())
virtual unsigned int order() const
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
unsigned int nSteps() const
double & A(unsigned int i, unsigned int j)
double & bHat(unsigned int i)
unsigned int order() const
double & b(unsigned int i)
std::vector< const AbsFunction * > _diffEqn