Chaste  Release::2018.1
RungeKuttaFehlbergIvpOdeSolver.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "RungeKuttaFehlbergIvpOdeSolver.hpp"
37 #include <cmath>
38 #include <cfloat>
39 #include <iostream>
40 #include "Exception.hpp"
41 
42 /*
43  * PROTECTED FUNCTIONS =========================================================
44  */
45 
47  AbstractOdeSystem* pOdeSystem,
48  std::vector<double>& rYValues,
49  std::vector<double>& rWorkingMemory,
50  double startTime,
51  double endTime,
52  double maxTimeStep,
53  double minTimeStep,
54  double tolerance,
55  bool outputSolution)
56 {
57  const unsigned number_of_variables = pOdeSystem->GetNumberOfStateVariables();
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);
70 
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;
76 
77  if (outputSolution)
78  { // Write out ICs
79  rSolution.rGetTimes().push_back(current_time);
80  rSolution.rGetSolutions().push_back(rYValues);
81  }
82 
83  // should never get here if this bool has been set to true;
84  assert(!mStoppingEventOccurred);
85  while (!got_to_end)
86  {
87  //std::cout << "New timestep\n" << std::flush;
88  while (!accurate_enough)
89  {
90  accurate_enough = true; // assume it is OK until we check and find otherwise
91 
92  // Function that calls the appropriate one-step solver
93  CalculateNextYValue(pOdeSystem,
94  time_step,
95  current_time,
96  rYValues,
97  rWorkingMemory);
98 
99  // Find the maximum error in this vector
100  double max_error = -DBL_MAX;
101  for (unsigned i=0; i<number_of_variables; i++)
102  {
103  if (mError[i] > max_error)
104  {
105  max_error = mError[i];
106  }
107  }
108 
109  if (max_error > tolerance)
110  {// Reject the step-size and do it again.
111  accurate_enough = false;
112  //std::cout << "Approximation rejected\n" << std::flush;
113  }
114  else
115  {
116  // step forward the time since step has now been made
117  current_time = current_time + time_step;
118  //std::cout << "Approximation accepted with time step = "<< time_step << "\n" << std::flush;
119  //std::cout << "max_error = " << max_error << " tolerance = " << tolerance << "\n" << std::flush;
120  if (outputSolution)
121  { // Write out ICs
122  //std::cout << "In solver Time = " << current_time << " y = " << rWorkingMemory[0] << "\n" << std::flush;
123  rSolution.rGetTimes().push_back(current_time);
124  rSolution.rGetSolutions().push_back(rWorkingMemory);
125  number_of_time_steps++;
126  }
127  }
128 
129  // Set a new step size based on the accuracy here
130  AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
131  }
132 
133  // For the next timestep check the step doesn't go past the end...
134  if (current_time + time_step > endTime)
135  { // Allow a smaller timestep for the final step.
136  time_step = endTime - current_time;
137  }
138 
139  if (pOdeSystem->CalculateStoppingEvent(current_time,
140  rWorkingMemory) == true)
141  {
142  mStoppingTime = current_time;
143  mStoppingEventOccurred = true;
144  }
145 
146  if (mStoppingEventOccurred || current_time>=endTime)
147  {
148  got_to_end = true;
149  }
150 
151  // Approximation accepted - put it in rYValues
152  rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
153  accurate_enough = false; // for next loop.
154  //std::cout << "Finished Time Step\n" << std::flush;
155  }
156  rSolution.SetNumberOfTimeSteps(number_of_time_steps);
157 }
158 
160  double timeStep,
161  double time,
162  std::vector<double>& rCurrentYValues,
163  std::vector<double>& rNextYValues)
164 {
165  const unsigned num_equations = pAbstractOdeSystem->GetNumberOfStateVariables();
166 
167 
168  std::vector<double>& dy = rNextYValues; // re-use memory (not that it makes much difference here!)
169 
170  pAbstractOdeSystem->EvaluateYDerivatives(time, rCurrentYValues, dy);
171 
172  for (unsigned i=0; i<num_equations; i++)
173  {
174  mk1[i] = timeStep*dy[i];
175  myk2[i] = rCurrentYValues[i] + 0.25*mk1[i];
176  }
177 
178  pAbstractOdeSystem->EvaluateYDerivatives(time + 0.25*timeStep, myk2, dy);
179  for (unsigned i=0; i<num_equations; i++)
180  {
181  mk2[i] = timeStep*dy[i];
182  myk3[i] = rCurrentYValues[i] + 0.09375*mk1[i] + 0.28125*mk2[i];
183  }
184 
185  pAbstractOdeSystem->EvaluateYDerivatives(time + 0.375*timeStep, myk3, dy);
186  for (unsigned i=0; i<num_equations; i++)
187  {
188  mk3[i] = timeStep*dy[i];
189  myk4[i] = rCurrentYValues[i] + m1932o2197*mk1[i] - m7200o2197*mk2[i]
190  + m7296o2197*mk3[i];
191  }
192 
193  pAbstractOdeSystem->EvaluateYDerivatives(time+m12o13*timeStep, myk4, dy);
194  for (unsigned i=0; i<num_equations; i++)
195  {
196  mk4[i] = timeStep*dy[i];
197  myk5[i] = rCurrentYValues[i] + m439o216*mk1[i] - 8*mk2[i]
198  + m3680o513*mk3[i]- m845o4104*mk4[i];
199  }
200 
201  pAbstractOdeSystem->EvaluateYDerivatives(time+timeStep, myk5, dy);
202  for (unsigned i=0; i<num_equations; i++)
203  {
204  mk5[i] = timeStep*dy[i];
205  myk6[i] = rCurrentYValues[i] - m8o27*mk1[i] + 2*mk2[i] - m3544o2565*mk3[i]
206  + m1859o4104*mk4[i] - 0.275*mk5[i];
207  }
208 
209  pAbstractOdeSystem->EvaluateYDerivatives(time+0.5*timeStep, myk6, dy);
210  for (unsigned i=0; i<num_equations; i++)
211  {
212  mk6[i] = timeStep*dy[i];
213  mError[i] = (1/timeStep)*fabs(m1o360*mk1[i] - m128o4275*mk3[i]
214  - m2197o75240*mk4[i] + 0.02*mk5[i]+ m2o55*mk6[i]);
215  rNextYValues[i] = rCurrentYValues[i] + m25o216*mk1[i] + m1408o2565*mk3[i]
216  + m2197o4104*mk4[i] - 0.2*mk5[i];
217  }
218 }
219 
221  const double& rError,
222  const double& rTolerance,
223  const double& rMaxTimeStep,
224  const double& rMinTimeStep)
225 {
226  // Work out scaling factor delta for the step size
227  double delta = pow(rTolerance/(2.0*rError), 0.25);
228 
229  // Maximum adjustment is *0.1 or *4
230  if (delta <= 0.1)
231  {
232  rCurrentStepSize *= 0.1;
233  }
234  else if (delta >= 4.0)
235  {
236  rCurrentStepSize *= 4.0;
237  }
238  else
239  {
240  rCurrentStepSize *= delta;
241  }
242 
243  if (rCurrentStepSize > rMaxTimeStep)
244  {
245  rCurrentStepSize = rMaxTimeStep;
246  }
247 
248  if (rCurrentStepSize < rMinTimeStep)
249  {
250  std::cout << "rCurrentStepSize = " << rCurrentStepSize << "\n" << std::flush;
251  std::cout << "rMinTimeStep = " << rMinTimeStep << "\n" << std::flush;
252 
253  EXCEPTION("RKF45 Solver: Ode needs a smaller timestep than the set minimum\n");
254  }
255 
256 }
257 
258 
259 /*
260  * PUBLIC FUNCTIONS=============================================================
261  */
262 
264  : m1932o2197(1932.0/2197.0), // you need the .0 s - caused me no end of trouble.
265  m7200o2197(7200.0/2197.0),
266  m7296o2197(7296.0/2197.0),
267  m12o13(12.0/13.0),
268  m439o216(439.0/216.0),
269  m3680o513(3680.0/513.0),
270  m845o4104(845.0/4104.0),
271  m8o27(8.0/27.0),
272  m3544o2565(3544.0/2565.0),
273  m1859o4104(1859.0/4104.0),
274  m1o360(1.0/360.0),
275  m128o4275(128.0/4275.0),
276  m2197o75240(2197.0/75240.0),
277  m2o55(2.0/55.0),
278  m25o216(25.0/216.0),
279  m1408o2565(1408.0/2565.0),
280  m2197o4104(2197.0/4104.0)
281 {
282 }
283 
285  std::vector<double>& rYValues,
286  double startTime,
287  double endTime,
288  double timeStep,
289  double tolerance)
290 {
291  assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
292  assert(endTime > startTime);
293  assert(timeStep > 0.0);
294 
295  mStoppingEventOccurred = false;
296  if (pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
297  {
298  EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
299  }
300  // Allocate working memory
301  std::vector<double> working_memory(rYValues.size());
302  // And solve...
303  OdeSolution solutions;
304  //solutions.SetNumberOfTimeSteps((unsigned)(10.0*(startTime-endTime)/timeStep));
305  bool return_solution = true;
306  InternalSolve(solutions, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-5, tolerance, return_solution);
307  return solutions;
308 }
309 
311  std::vector<double>& rYValues,
312  double startTime,
313  double endTime,
314  double timeStep)
315 {
316  assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
317  assert(endTime > startTime);
318  assert(timeStep > 0.0);
319 
320  mStoppingEventOccurred = false;
321  if (pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
322  {
323  EXCEPTION("(Solve without sampling) Stopping event is true for initial condition");
324  }
325  // Allocate working memory
326  std::vector<double> working_memory(rYValues.size());
327  // And solve...
328  OdeSolution not_required_solution;
329  bool return_solution = false;
330  InternalSolve(not_required_solution, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-4, 1e-5, return_solution);
331 }
332 
333 
334 // Serialization for Boost >= 1.36
std::vector< std::vector< double > > & rGetSolutions()
#define EXCEPTION(message)
Definition: Exception.hpp:143
void CalculateNextYValue(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rNextYValues)
virtual void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)=0
OdeSolution Solve(AbstractOdeSystem *pAbstractOdeSystem, std::vector< double > &rYValues, double startTime, double endTime, double timeStep, double ignoredSamplingTime)
std::vector< double > & rGetTimes()
void AdjustStepSize(double &rCurrentStepSize, const double &rError, const double &rTolerance, const double &rMaxTimeStep, const double &rMinTimeStep)
#define CHASTE_CLASS_EXPORT(T)
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)
void SetNumberOfTimeSteps(unsigned numTimeSteps)
Definition: OdeSolution.cpp:58
virtual bool CalculateStoppingEvent(double time, const std::vector< double > &rY)