48 std::vector<double>& rYValues,
49 std::vector<double>& rWorkingMemory,
58 mError.resize(number_of_variables);
59 mk1.resize(number_of_variables);
60 mk2.resize(number_of_variables);
61 mk3.resize(number_of_variables);
62 mk4.resize(number_of_variables);
63 mk5.resize(number_of_variables);
64 mk6.resize(number_of_variables);
65 myk2.resize(number_of_variables);
66 myk3.resize(number_of_variables);
67 myk4.resize(number_of_variables);
68 myk5.resize(number_of_variables);
69 myk6.resize(number_of_variables);
71 double current_time = startTime;
72 double time_step = maxTimeStep;
73 bool got_to_end =
false;
74 bool accurate_enough =
false;
75 unsigned number_of_time_steps = 0;
79 rSolution.
rGetTimes().push_back(current_time);
88 while (!accurate_enough)
90 accurate_enough =
true;
100 double max_error = -DBL_MAX;
101 for (
unsigned i=0; i<number_of_variables; i++)
103 if (
mError[i] > max_error)
109 if (max_error > tolerance)
111 accurate_enough =
false;
117 current_time = current_time + time_step;
123 rSolution.
rGetTimes().push_back(current_time);
125 number_of_time_steps++;
130 AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
134 if (current_time + time_step > endTime)
136 time_step = endTime - current_time;
140 rWorkingMemory) ==
true)
152 rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
153 accurate_enough =
false;
162 std::vector<double>& rCurrentYValues,
163 std::vector<double>& rNextYValues)
168 std::vector<double>& dy = rNextYValues;
172 for (
unsigned i=0; i<num_equations; i++)
174 mk1[i] = timeStep*dy[i];
175 myk2[i] = rCurrentYValues[i] + 0.25*
mk1[i];
179 for (
unsigned i=0; i<num_equations; i++)
181 mk2[i] = timeStep*dy[i];
182 myk3[i] = rCurrentYValues[i] + 0.09375*
mk1[i] + 0.28125*
mk2[i];
186 for (
unsigned i=0; i<num_equations; i++)
188 mk3[i] = timeStep*dy[i];
194 for (
unsigned i=0; i<num_equations; i++)
196 mk4[i] = timeStep*dy[i];
202 for (
unsigned i=0; i<num_equations; i++)
204 mk5[i] = timeStep*dy[i];
210 for (
unsigned i=0; i<num_equations; i++)
212 mk6[i] = timeStep*dy[i];
221 const double& rError,
222 const double& rTolerance,
223 const double& rMaxTimeStep,
224 const double& rMinTimeStep)
227 double delta = pow(rTolerance/(2.0*rError), 0.25);
232 rCurrentStepSize *= 0.1;
234 else if (delta >= 4.0)
236 rCurrentStepSize *= 4.0;
240 rCurrentStepSize *= delta;
243 if (rCurrentStepSize > rMaxTimeStep)
245 rCurrentStepSize = rMaxTimeStep;
248 if (rCurrentStepSize < rMinTimeStep)
250 std::cout <<
"rCurrentStepSize = " << rCurrentStepSize <<
"\n" << std::flush;
251 std::cout <<
"rMinTimeStep = " << rMinTimeStep <<
"\n" << std::flush;
253 EXCEPTION(
"RKF45 Solver: Ode needs a smaller timestep than the set minimum\n");
264 : m1932o2197(1932.0/2197.0),
265 m7200o2197(7200.0/2197.0),
266 m7296o2197(7296.0/2197.0),
268 m439o216(439.0/216.0),
269 m3680o513(3680.0/513.0),
270 m845o4104(845.0/4104.0),
272 m3544o2565(3544.0/2565.0),
273 m1859o4104(1859.0/4104.0),
275 m128o4275(128.0/4275.0),
276 m2197o75240(2197.0/75240.0),
279 m1408o2565(1408.0/2565.0),
280 m2197o4104(2197.0/4104.0)
285 std::vector<double>& rYValues,
292 assert(endTime > startTime);
293 assert(timeStep > 0.0);
298 EXCEPTION(
"(Solve with sampling) Stopping event is true for initial condition");
301 std::vector<double> working_memory(rYValues.size());
305 bool return_solution =
true;
306 InternalSolve(solutions, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-5, tolerance, return_solution);
311 std::vector<double>& rYValues,
317 assert(endTime > startTime);
318 assert(timeStep > 0.0);
323 EXCEPTION(
"(Solve without sampling) Stopping event is true for initial condition");
326 std::vector<double> working_memory(rYValues.size());
329 bool return_solution =
false;
330 InternalSolve(not_required_solution, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-4, 1e-5, return_solution);
void InternalSolve(OdeSolution &rSolution, AbstractOdeSystem *pAbstractOdeSystem, std::vector< double > &rCurrentYValues, std::vector< double > &rWorkingMemory, double startTime, double endTime, double maxTimeStep, double minTimeStep, double tolerance, bool outputSolution)