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>
4namespace Genfun {
5
6
7 StepDoublingRKStepper::StepDoublingRKStepper(const ButcherTableau & xtableau):tableau(xtableau) {
8 }
9
11 }
12
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
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}
unsigned int nSteps() const
double & A(unsigned int i, unsigned int j)
double & b(unsigned int i)
unsigned int order() const
std::vector< const AbsFunction * > _diffEqn
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, std::vector< double > &errors) const
virtual unsigned int order() const
void doStep(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &s, RKIntegrator::RKData::Data &d) const
virtual StepDoublingRKStepper * clone() const
StepDoublingRKStepper(const ButcherTableau &tableau)