CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 #include "CLHEP/GenericFunctions/EmbeddedRKStepper.hh" 00002 #include "CLHEP/GenericFunctions/ExtendedButcherTableau.hh" 00003 #include <stdexcept> 00004 namespace Genfun { 00005 00006 00007 EmbeddedRKStepper::EmbeddedRKStepper(const ExtendedButcherTableau & mtableau): 00008 tableau(mtableau){ 00009 } 00010 00011 EmbeddedRKStepper::~EmbeddedRKStepper() { 00012 } 00013 00014 void EmbeddedRKStepper::step (const RKIntegrator::RKData * data, 00015 const RKIntegrator::RKData::Data & s, 00016 RKIntegrator::RKData::Data & d, 00017 std::vector<double> & errors) const { 00018 00019 // First step: 00020 double h = d.time - s.time; 00021 if (h<=0) throw std::runtime_error ("Runtime error in RKIntegrator (zero or negative stepsize)"); 00022 unsigned int nvar = s.variable.size(); 00023 00024 // Compute all of the k's..: 00025 // 00026 std::vector<std::vector<double> >k(tableau.nSteps()); 00027 for (unsigned int i=0;i<tableau.nSteps();i++) { 00028 k[i].resize(nvar,0); 00029 Argument arg(nvar); 00030 for (unsigned int v=0;v<nvar;v++) arg[v]=s.variable[v]; 00031 for (unsigned int j=0;j<i;j++) { 00032 for (unsigned int v=0;v<nvar;v++) arg[v] += h*tableau.A(i,j)*k[j][v]; 00033 } 00034 for (unsigned int v=0;v<nvar;v++) k[i][v]=(*data->_diffEqn[v])(arg); 00035 } 00036 // 00037 // Final result. 00038 // 00039 for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] = 0; 00040 for (unsigned int i=0;i<tableau.nSteps();i++) { 00041 for (unsigned int v=0;v<nvar;v++) d.firstDerivative[v] += tableau.b(i)*k[i][v]; 00042 } 00043 for (unsigned int v=0;v<nvar;v++) d.variable[v] =s.variable[v]+h*d.firstDerivative[v]; 00044 // 00045 // 00046 // 00047 errors.resize(nvar); 00048 for (unsigned int v=0;v<nvar;v++) errors[v] = 0; 00049 for (unsigned int i=0;i<tableau.nSteps();i++) { 00050 for (unsigned int v=0;v<nvar;v++) errors[v] += (h*(tableau.bHat(i)-tableau.b(i))*k[i][v]); 00051 } 00052 return; 00053 } 00054 00055 EmbeddedRKStepper *EmbeddedRKStepper::clone() const { 00056 return new EmbeddedRKStepper(*this); 00057 } 00058 00059 unsigned int EmbeddedRKStepper::order() const { 00060 return tableau.order(); 00061 } 00062 }