00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "RungeKuttaFehlbergIvpOdeSolver.hpp"
00030 #include <cmath>
00031 #include <cfloat>
00032 #include <iostream>
00033 #include "Exception.hpp"
00034
00035 const double smidge = 1e-10;
00036
00037
00038
00039
00040
00041 void RungeKuttaFehlbergIvpOdeSolver::InternalSolve(OdeSolution& rSolution,
00042 AbstractOdeSystem* pOdeSystem,
00043 std::vector<double>& rYValues,
00044 std::vector<double>& rWorkingMemory,
00045 double startTime,
00046 double endTime,
00047 double maxTimeStep,
00048 double minTimeStep,
00049 double tolerance,
00050 bool outputSolution)
00051 {
00052 const unsigned number_of_variables = pOdeSystem->GetNumberOfStateVariables();
00053 mError.reserve(number_of_variables);
00054 mk1.reserve(number_of_variables);
00055 mk2.reserve(number_of_variables);
00056 mk3.reserve(number_of_variables);
00057 mk4.reserve(number_of_variables);
00058 mk5.reserve(number_of_variables);
00059 mk6.reserve(number_of_variables);
00060 myk2.reserve(number_of_variables);
00061 myk3.reserve(number_of_variables);
00062 myk4.reserve(number_of_variables);
00063 myk5.reserve(number_of_variables);
00064 myk6.reserve(number_of_variables);
00065
00066 double current_time = startTime;
00067 double time_step = maxTimeStep;
00068 bool got_to_end = false;
00069 bool accurate_enough = false;
00070 unsigned number_of_time_steps = 0;
00071
00072 if (outputSolution)
00073 {
00074 rSolution.rGetTimes().push_back(current_time);
00075 rSolution.rGetSolutions().push_back(rYValues);
00076 }
00077
00078
00079 assert(!mStoppingEventOccurred);
00080 while (!got_to_end)
00081 {
00082
00083 while (!accurate_enough)
00084 {
00085 accurate_enough = true;
00086
00087
00088 CalculateNextYValue(pOdeSystem,
00089 time_step,
00090 current_time,
00091 rYValues,
00092 rWorkingMemory);
00093
00094
00095 double max_error = -DBL_MAX;
00096 for (unsigned i=0; i<number_of_variables; i++)
00097 {
00098 if (mError[i] > max_error)
00099 {
00100 max_error = mError[i];
00101 }
00102 }
00103
00104 if (max_error > tolerance)
00105 {
00106 accurate_enough = false;
00107
00108 }
00109 else
00110 {
00111
00112 current_time = current_time + time_step;
00113
00114
00115 if (outputSolution)
00116 {
00117
00118 rSolution.rGetTimes().push_back(current_time);
00119 rSolution.rGetSolutions().push_back(rWorkingMemory);
00120 number_of_time_steps++;
00121 }
00122 }
00123
00124
00125 AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
00126 }
00127
00128
00129 if (current_time + time_step > endTime)
00130 {
00131 time_step = endTime - current_time;
00132 }
00133
00134 if ( pOdeSystem->CalculateStoppingEvent(current_time,
00135 rWorkingMemory) == true )
00136 {
00137 mStoppingTime = current_time;
00138 mStoppingEventOccurred = true;
00139 }
00140
00141 if (mStoppingEventOccurred || current_time>=endTime)
00142 {
00143 got_to_end = true;
00144 }
00145
00146
00147 rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
00148 accurate_enough = false;
00149
00150 }
00151 rSolution.SetNumberOfTimeSteps(number_of_time_steps);
00152 }
00153
00154 void RungeKuttaFehlbergIvpOdeSolver::CalculateNextYValue(AbstractOdeSystem* pAbstractOdeSystem,
00155 double timeStep,
00156 double time,
00157 std::vector<double>& rCurrentYValues,
00158 std::vector<double>& rNextYValues)
00159 {
00160 const unsigned num_equations = pAbstractOdeSystem->GetNumberOfStateVariables();
00161
00162
00163 std::vector<double>& dy = rNextYValues;
00164
00165 pAbstractOdeSystem->EvaluateYDerivatives(time, rCurrentYValues, dy);
00166
00167 for (unsigned i=0; i<num_equations; i++)
00168 {
00169 mk1[i] = timeStep*dy[i];
00170 myk2[i] = rCurrentYValues[i] + 0.25*mk1[i];
00171 }
00172
00173 pAbstractOdeSystem->EvaluateYDerivatives(time + 0.25*timeStep, myk2, dy);
00174 for (unsigned i=0; i<num_equations; i++)
00175 {
00176 mk2[i] = timeStep*dy[i];
00177 myk3[i] = rCurrentYValues[i] + 0.09375*mk1[i] + 0.28125*mk2[i];
00178 }
00179
00180 pAbstractOdeSystem->EvaluateYDerivatives(time + 0.375*timeStep, myk3, dy);
00181 for (unsigned i=0; i<num_equations; i++)
00182 {
00183 mk3[i] = timeStep*dy[i];
00184 myk4[i] = rCurrentYValues[i] + m1932o2197*mk1[i] - m7200o2197*mk2[i]
00185 + m7296o2197*mk3[i];
00186 }
00187
00188 pAbstractOdeSystem->EvaluateYDerivatives(time+m12o13*timeStep, myk4, dy);
00189 for (unsigned i=0; i<num_equations; i++)
00190 {
00191 mk4[i] = timeStep*dy[i];
00192 myk5[i] = rCurrentYValues[i] + m439o216*mk1[i] - 8*mk2[i]
00193 + m3680o513*mk3[i]- m845o4104*mk4[i];
00194 }
00195
00196 pAbstractOdeSystem->EvaluateYDerivatives(time+timeStep, myk5, dy);
00197 for (unsigned i=0; i<num_equations; i++)
00198 {
00199 mk5[i] = timeStep*dy[i];
00200 myk6[i] = rCurrentYValues[i] - m8o27*mk1[i] + 2*mk2[i] - m3544o2565*mk3[i]
00201 + m1859o4104*mk4[i] - 0.275*mk5[i];
00202 }
00203
00204 pAbstractOdeSystem->EvaluateYDerivatives(time+0.5*timeStep, myk6, dy);
00205 for (unsigned i=0; i<num_equations; i++)
00206 {
00207 mk6[i] = timeStep*dy[i];
00208 mError[i] = (1/timeStep)*fabs(m1o360*mk1[i] - m128o4275*mk3[i]
00209 - m2197o75240*mk4[i] + 0.02*mk5[i]+ m2o55*mk6[i]);
00210 rNextYValues[i] = rCurrentYValues[i] + m25o216*mk1[i] + m1408o2565*mk3[i]
00211 + m2197o4104*mk4[i] - 0.2*mk5[i];
00212 }
00213 }
00214
00215 void RungeKuttaFehlbergIvpOdeSolver::AdjustStepSize(double& rCurrentStepSize,
00216 const double& rError,
00217 const double& rTolerance,
00218 const double& rMaxTimeStep,
00219 const double& rMinTimeStep)
00220 {
00221
00222 double delta = pow(rTolerance/(2.0*rError), 0.25);
00223
00224
00225 if (delta <= 0.1)
00226 {
00227 rCurrentStepSize *= 0.1;
00228 }
00229 else if (delta >= 4.0)
00230 {
00231 rCurrentStepSize *= 4.0;
00232 }
00233 else
00234 {
00235 rCurrentStepSize *= delta;
00236 }
00237
00238 if (rCurrentStepSize > rMaxTimeStep)
00239 {
00240 rCurrentStepSize = rMaxTimeStep;
00241 }
00242
00243 if (rCurrentStepSize < rMinTimeStep)
00244 {
00245 std::cout << "rCurrentStepSize = " << rCurrentStepSize << "\n" << std::flush;
00246 std::cout << "rMinTimeStep = " << rMinTimeStep << "\n" << std::flush;
00247
00248 EXCEPTION("RKF45 Solver: Ode needs a smaller timestep than the set minimum\n");
00249 }
00250
00251 }
00252
00253
00254
00255
00256
00257
00258 RungeKuttaFehlbergIvpOdeSolver::RungeKuttaFehlbergIvpOdeSolver()
00259 : m1932o2197(1932.0/2197.0),
00260 m7200o2197(7200.0/2197.0),
00261 m7296o2197(7296.0/2197.0),
00262 m12o13(12.0/13.0),
00263 m439o216(439.0/216.0),
00264 m3680o513(3680.0/513.0),
00265 m845o4104(845.0/4104.0),
00266 m8o27(8.0/27.0),
00267 m3544o2565(3544.0/2565.0),
00268 m1859o4104(1859.0/4104.0),
00269 m1o360(1.0/360.0),
00270 m128o4275(128.0/4275.0),
00271 m2197o75240(2197.0/75240.0),
00272 m2o55(2.0/55.0),
00273 m25o216(25.0/216.0),
00274 m1408o2565(1408.0/2565.0),
00275 m2197o4104(2197.0/4104.0)
00276 {
00277 }
00278
00279 OdeSolution RungeKuttaFehlbergIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00280 std::vector<double>& rYValues,
00281 double startTime,
00282 double endTime,
00283 double timeStep,
00284 double tolerance)
00285 {
00286 assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00287 assert(endTime > startTime);
00288 assert(timeStep > 0.0);
00289
00290 mStoppingEventOccurred = false;
00291 if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00292 {
00293 EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
00294 }
00295
00296 std::vector<double> working_memory(rYValues.size());
00297
00298 OdeSolution solutions;
00299
00300 bool return_solution = true;
00301 InternalSolve(solutions, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-5, tolerance, return_solution);
00302 return solutions;
00303 }
00304
00305 void RungeKuttaFehlbergIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00306 std::vector<double>& rYValues,
00307 double startTime,
00308 double endTime,
00309 double timeStep)
00310 {
00311 assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00312 assert(endTime > startTime);
00313 assert(timeStep > 0.0);
00314
00315 mStoppingEventOccurred = false;
00316 if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00317 {
00318 EXCEPTION("(Solve without sampling) Stopping event is true for initial condition");
00319 }
00320
00321 std::vector<double> working_memory(rYValues.size());
00322
00323 OdeSolution not_required_solution;
00324 bool return_solution = false;
00325 InternalSolve(not_required_solution, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-4, 1e-5, return_solution);
00326 }
00327
00328
00329
00330 #include "SerializationExportWrapperForCpp.hpp"
00331 CHASTE_CLASS_EXPORT(RungeKuttaFehlbergIvpOdeSolver)