RungeKuttaFehlbergIvpOdeSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
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  * PROTECTED FUNCTIONS =========================================================
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     {   // Write out ICs
00074         rSolution.rGetTimes().push_back(current_time);
00075         rSolution.rGetSolutions().push_back(rYValues);
00076     }
00077 
00078     // should never get here if this bool has been set to true;
00079     assert(!mStoppingEventOccurred);
00080     while (!got_to_end)
00081     {
00082         //std::cout << "New timestep\n" << std::flush;
00083         while (!accurate_enough)
00084         {
00085             accurate_enough = true; // assume it is OK until we check and find otherwise
00086 
00087             // Function that calls the appropriate one-step solver
00088             CalculateNextYValue(pOdeSystem,
00089                                 time_step,
00090                                 current_time,
00091                                 rYValues,
00092                                 rWorkingMemory);
00093 
00094             // Find the maximum error in this vector
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             {// Reject the step-size and do it again.
00106                 accurate_enough = false;
00107                 //std::cout << "Approximation rejected\n" << std::flush;
00108             }
00109             else
00110             {
00111                 // step forward the time since step has now been made
00112                 current_time = current_time + time_step;
00113                 //std::cout << "Approximation accepted with time step = "<< time_step << "\n" << std::flush;
00114                 //std::cout << "max_error = " << max_error << " tolerance = " << tolerance << "\n" << std::flush;
00115                 if (outputSolution)
00116                 {   // Write out ICs
00117                     //std::cout << "In solver Time = " << current_time << " y = " << rWorkingMemory[0] << "\n" << std::flush;
00118                     rSolution.rGetTimes().push_back(current_time);
00119                     rSolution.rGetSolutions().push_back(rWorkingMemory);
00120                     number_of_time_steps++;
00121                 }
00122             }
00123 
00124             // Set a new step size based on the accuracy here
00125             AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
00126         }
00127 
00128         // For the next timestep check the step doesn't go past the end...
00129         if (current_time + time_step > endTime)
00130         {   // Allow a smaller timestep for the final step.
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         // Approximation accepted - put it in rYValues
00147         rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
00148         accurate_enough = false; // for next loop.
00149         //std::cout << "Finished Time Step\n" << std::flush;
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; // re-use memory (not that it makes much difference here!)
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     // Work out scaling factor delta for the step size
00222     double delta = pow(rTolerance/(2.0*rError), 0.25);
00223 
00224     // Maximum adjustment is *0.1 or *4
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  * PUBLIC FUNCTIONS=============================================================
00256  */
00257 
00258 RungeKuttaFehlbergIvpOdeSolver::RungeKuttaFehlbergIvpOdeSolver()
00259     : m1932o2197(1932.0/2197.0), // you need the .0 s - caused me no end of trouble.
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     // Allocate working memory
00296     std::vector<double> working_memory(rYValues.size());
00297     // And solve...
00298     OdeSolution solutions;
00299     //solutions.SetNumberOfTimeSteps((unsigned)(10.0*(startTime-endTime)/timeStep));
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     // Allocate working memory
00321     std::vector<double> working_memory(rYValues.size());
00322     // And solve...
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 // Serialization for Boost >= 1.36
00330 #include "SerializationExportWrapperForCpp.hpp"
00331 CHASTE_CLASS_EXPORT(RungeKuttaFehlbergIvpOdeSolver)

Generated on Mon Apr 18 11:35:36 2011 for Chaste by  doxygen 1.5.5