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