42 if (_data->
_fx.empty()) {
46 for (
size_t f=0;
f<nvar;
f++) {
56 std::set<RKData::Data>::iterator s=_data->
_fx.find(dt);
57 if (s!=_data->
_fx.end()) {
60 return (*s).variable[_index];
63 s=_data->
_fx.lower_bound(dt);
65 if (!(s!=_data->
_fx.begin()))
throw std::runtime_error(
"Runtime error in RKIntegrator");
70 _data->
_fx.insert(s,dt);
80 RKIntegrator::RKData::~RKData() {
81 for (
size_t i=0;i<_startingValParameter.size();i++)
delete _startingValParameter[i];
82 for (
size_t i=0;i<_controlParameter.size();i++)
delete _controlParameter[i];
83 for (
size_t i=0;i<_diffEqn.size(); i++)
delete _diffEqn[i];
94 for (
size_t i=0;i<_fcn.size();i++)
delete _fcn[i];
98 const std::string &variableName,
99 double defStartingValue,
101 double defValueMax) {
102 Parameter *par =
new Parameter(variableName, defStartingValue, defValueMin, defValueMax);
106 _fcn.push_back(
new RKFunction(_data,_fcn.size()));
115 double defStartingValue,
116 double startingValueMin,
117 double startingValueMax) {
119 Parameter *par =
new Parameter(variableName, defStartingValue, startingValueMin, startingValueMax);
136 unsigned int size =
_diffEqn.size();
137 for (
size_t i=0;i<size;i++) {
138 if (!(
_diffEqn[i]->dimensionality()==size))
throw std::runtime_error(
"Runtime error in RKIntegrator");
150 for (
size_t p=0;p<_startingValParameter.size();p++) {
151 if (_startingValParameter[p]->getValue()!=_startingValParameterCache[p]) {
152 _startingValParameterCache[p]=_startingValParameter[p]->getValue();
160 for (
size_t p=0;p<_controlParameter.size();p++) {
161 if (_controlParameter[p]->getValue()!=_controlParameterCache[p]) {
162 _controlParameterCache[p]=_controlParameter[p]->getValue();
170 _fx.erase(_fx.begin(),_fx.end());
182 const double eps = 1.0E-6;
183 const double SAFETY = 0.9;
184 const double PSHRNK = -0.25;
185 const double PGROW = -0.20;
186 const double ERRCON = -1.89E-4;
187 const double TINY = 1.0E-30;
201 double h = Tmp1.
time - Tmp0.time;
203 std::vector<double> errors;
205 rkck(Tmp0, Tmp1, errors);
206 for (
size_t e=0;e<errors.size();e++) {
207 errors[e] = fabs(errors[e]) / (fabs(Tmp0.variable[e]) + fabs(Tmp0.firstDerivative[e]*h) + TINY);
209 double emax = (*std::max_element(errors.begin(),errors.end()))/eps;
211 h = std::max(SAFETY*h*pow(emax,PSHRNK),0.1*h);
212 if (!(((
float) Tmp0.time+h - (
float) Tmp0.time) > 0) ) {
213 std::cerr <<
"Warning, RK Integrator step underflow" << std::endl;
215 Tmp1.time = Tmp0.time+h;
221 hnext = SAFETY*h*pow(emax,PGROW);
232 Tmp1.time = std::min(Tmp0.time + hnext, d.
time);
247 void RKIntegrator::RKFunction::rkck(
const RKIntegrator::RKData::Data & s, RKIntegrator::RKData::Data & d, std::vector<double> & errors)
const {
249 #ifdef NONAUTONOMOUS_EQUATIONS
272 b64=44275.0/110592.0,
282 dc1=c1-2825.0/27648.0,
283 dc3=c3-18575.0/48384.0,
284 dc4=c4-13525.0/55296.0,
289 double h = d.
time - s.time;
290 if (h<=0)
throw std::runtime_error (
"Runtime error in RKIntegrator (zero or negative stepsize)");
291 unsigned int nv = s.variable.size();
292 Argument arg(nv), arg0(nv), d1(nv),d2(nv), d3(nv), d4(nv), d5(nv), d6(nv);
295 for (
size_t v=0;v<nv;v++) { arg0[v]=s.variable[v];}
298 for (
size_t v=0;v<nv;v++) {d1[v]=(*_data->
_diffEqn[v])(arg0);}
299 for (
size_t v=0;v<nv;v++) {s.firstDerivative[v]=d1[v];}
303 for (
size_t v=0;v<nv;v++) { d1[v] = s.firstDerivative[v];}
307 for (
size_t v=0;v<nv;v++) { arg[v] = arg0[v] + b21*h*d1[v];}
309 for (
size_t v=0;v<nv;v++) { d2[v] = (*_data->
_diffEqn[v])(arg);}
310 for (
size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b31*d1[v]+b32*d2[v]);}
313 for (
size_t v=0;v<nv;v++) { d3[v] = (*_data->
_diffEqn[v])(arg);}
314 for (
size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b41*d1[v]+b42*d2[v]+b43*d3[v]);}
316 for (
size_t v=0;v<nv;v++) { d4[v] = (*_data->
_diffEqn[v])(arg);}
317 for (
size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b51*d1[v]+b52*d2[v]+b53*d3[v] + b54*d4[v]);}
319 for (
size_t v=0;v<nv;v++) { d5[v] = (*_data->
_diffEqn[v])(arg);}
320 for (
size_t v=0;v<nv;v++) { arg[v] = arg0[v] + h*(b61*d1[v]+b62*d2[v]+b63*d3[v] + b64*d4[v] + b65*d5[v]);}
322 for (
size_t v=0;v<nv;v++) { d6[v] = (*_data->
_diffEqn[v])(arg);}
324 for (
size_t v=0;v<nv;v++) { d.variable[v] = arg0[v] + h*(c1*d1[v]+c3*d3[v]+c4*d4[v]+c6*d6[v]);}
325 errors.erase(errors.begin(),errors.end());
327 for (
size_t v=0;v<nv;v++) { errors.push_back(h*(dc1*d1[v]+dc3*d3[v]+dc4*d4[v]+dc5*d5[v]+dc6*d6[v]));}