17 std::vector<double> & errors)
const {
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();
26 std::vector<std::vector<double> >k(tableau.
nSteps());
27 for (
unsigned int i=0;i<tableau.
nSteps();i++) {
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];
34 for (
unsigned int v=0;v<nvar;v++) k[i][v]=(*data->
_diffEqn[v])(arg);
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];
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]);
60 return tableau.
order();
EmbeddedRKStepper(const ExtendedButcherTableau &tableau=CashKarpXtTableau())
virtual ~EmbeddedRKStepper()
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
std::vector< double > firstDerivative
std::vector< double > variable