Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
RungeKuttaFehlbergIvpOdeSolver.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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;
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;
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
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
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
#define EXCEPTION(message)
#define CHASTE_CLASS_EXPORT(T)
virtual void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)=0
virtual bool CalculateStoppingEvent(double time, const std::vector< double > &rY)
void SetNumberOfTimeSteps(unsigned numTimeSteps)
std::vector< std::vector< double > > & rGetSolutions()
std::vector< double > & rGetTimes()
void AdjustStepSize(double &rCurrentStepSize, const double &rError, const double &rTolerance, const double &rMaxTimeStep, const double &rMinTimeStep)
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)
OdeSolution Solve(AbstractOdeSystem *pAbstractOdeSystem, std::vector< double > &rYValues, double startTime, double endTime, double timeStep, double ignoredSamplingTime)
void CalculateNextYValue(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rNextYValues)