RungeKuttaFehlbergIvpOdeSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 
00034 const double smidge = 1e-10;
00035 
00036 /*
00037  * PROTECTED FUNCTIONS =========================================================
00038  */
00039 
00040 void RungeKuttaFehlbergIvpOdeSolver::InternalSolve(OdeSolution& rSolution,
00041                                                 AbstractOdeSystem* pOdeSystem,
00042                                                 std::vector<double>& rYValues,
00043                                                 std::vector<double>& rWorkingMemory,
00044                                                 double startTime,
00045                                                 double endTime,
00046                                                 double maxTimeStep,
00047                                                 double minTimeStep,
00048                                                 double tolerance,
00049                                                 bool outputSolution)
00050 {
00051     const unsigned number_of_variables = pOdeSystem->GetNumberOfStateVariables();
00052     mError.reserve(number_of_variables);
00053     mk1.reserve(number_of_variables);
00054     mk2.reserve(number_of_variables);
00055     mk3.reserve(number_of_variables);
00056     mk4.reserve(number_of_variables);
00057     mk5.reserve(number_of_variables);
00058     mk6.reserve(number_of_variables);
00059     myk2.reserve(number_of_variables);
00060     myk3.reserve(number_of_variables);
00061     myk4.reserve(number_of_variables);
00062     myk5.reserve(number_of_variables);
00063     myk6.reserve(number_of_variables);
00064 
00065     double current_time = startTime;
00066     double time_step = maxTimeStep;
00067     bool got_to_end = false;
00068     bool accurate_enough = false;
00069     unsigned number_of_time_steps = 0;
00070 
00071     if (outputSolution)
00072     {   // Write out ICs
00073         rSolution.rGetTimes().push_back(current_time);
00074         rSolution.rGetSolutions().push_back(rYValues);
00075     }
00076 
00077     // should never get here if this bool has been set to true;
00078     assert(!mStoppingEventOccurred);
00079     while (!got_to_end)
00080     {
00081         //std::cout << "New timestep\n" << std::flush;
00082         while (!accurate_enough)
00083         {
00084             accurate_enough = true; // assume it is OK until we check and find otherwise
00085 
00086             // Function that calls the appropriate one-step solver
00087             CalculateNextYValue(pOdeSystem,
00088                                 time_step,
00089                                 current_time,
00090                                 rYValues,
00091                                 rWorkingMemory);
00092 
00093             // Find the maximum error in this vector
00094             double max_error = -DBL_MAX;
00095             for (unsigned i=0; i<number_of_variables; i++)
00096             {
00097                 if (mError[i] > max_error)
00098                 {
00099                     max_error = mError[i];
00100                 }
00101             }
00102 
00103             if (max_error > tolerance)
00104             {// Reject the step-size and do it again.
00105                 accurate_enough = false;
00106                 //std::cout << "Approximation rejected\n" << std::flush;
00107             }
00108             else
00109             {
00110                 // step forward the time since step has now been made
00111                 current_time = current_time + time_step;
00112                 //std::cout << "Approximation accepted with time step = "<< time_step << "\n" << std::flush;
00113                 //std::cout << "max_error = " << max_error << " tolerance = " << tolerance << "\n" << std::flush;
00114                 if (outputSolution)
00115                 {   // Write out ICs
00116                     //std::cout << "In solver Time = " << current_time << " y = " << rWorkingMemory[0] << "\n" << std::flush;
00117                     rSolution.rGetTimes().push_back(current_time);
00118                     rSolution.rGetSolutions().push_back(rWorkingMemory);
00119                     number_of_time_steps++;
00120                 }
00121             }
00122 
00123             // Set a new step size based on the accuracy here
00124             AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
00125         }
00126 
00127         // For the next timestep check the step doesn't go past the end...
00128         if (current_time + time_step > endTime)
00129         {   // Allow a smaller timestep for the final step.
00130             time_step = endTime - current_time;
00131         }
00132 
00133         if ( pOdeSystem->CalculateStoppingEvent(current_time,
00134                                                 rWorkingMemory) == true )
00135         {
00136             mStoppingTime = current_time;
00137             mStoppingEventOccurred = true;
00138         }
00139 
00140         if (mStoppingEventOccurred || current_time>=endTime)
00141         {
00142             got_to_end = true;
00143         }
00144 
00145         // Approximation accepted - put it in rYValues
00146         rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
00147         accurate_enough = false; // for next loop.
00148         //std::cout << "Finished Time Step\n" << std::flush;
00149     }
00150     rSolution.SetNumberOfTimeSteps(number_of_time_steps);
00151 }
00152 
00153 void RungeKuttaFehlbergIvpOdeSolver::CalculateNextYValue(AbstractOdeSystem* pAbstractOdeSystem,
00154                                                   double timeStep,
00155                                                   double time,
00156                                                   std::vector<double>& rCurrentYValues,
00157                                                   std::vector<double>& rNextYValues)
00158 {
00159     const unsigned num_equations = pAbstractOdeSystem->GetNumberOfStateVariables();
00160 
00161 
00162     std::vector<double>& dy = rNextYValues; // re-use memory (not that it makes much difference here!)
00163 
00164     pAbstractOdeSystem->EvaluateYDerivatives(time, rCurrentYValues, dy);
00165 
00166     for (unsigned i=0; i<num_equations; i++)
00167     {
00168         mk1[i] = timeStep*dy[i];
00169         myk2[i] = rCurrentYValues[i] + 0.25*mk1[i];
00170     }
00171 
00172     pAbstractOdeSystem->EvaluateYDerivatives(time + 0.25*timeStep, myk2, dy);
00173     for (unsigned i=0; i<num_equations; i++)
00174     {
00175         mk2[i] = timeStep*dy[i];
00176         myk3[i] = rCurrentYValues[i] + 0.09375*mk1[i] + 0.28125*mk2[i];
00177     }
00178 
00179     pAbstractOdeSystem->EvaluateYDerivatives(time + 0.375*timeStep, myk3, dy);
00180     for (unsigned i=0; i<num_equations; i++)
00181     {
00182         mk3[i] = timeStep*dy[i];
00183         myk4[i] = rCurrentYValues[i] + m1932o2197*mk1[i] - m7200o2197*mk2[i]
00184                     + m7296o2197*mk3[i];
00185     }
00186 
00187     pAbstractOdeSystem->EvaluateYDerivatives(time+m12o13*timeStep, myk4, dy);
00188     for (unsigned i=0; i<num_equations; i++)
00189     {
00190         mk4[i] = timeStep*dy[i];
00191         myk5[i] = rCurrentYValues[i] + m439o216*mk1[i] - 8*mk2[i]
00192                     + m3680o513*mk3[i]- m845o4104*mk4[i];
00193     }
00194 
00195     pAbstractOdeSystem->EvaluateYDerivatives(time+timeStep, myk5, dy);
00196     for (unsigned i=0; i<num_equations; i++)
00197     {
00198         mk5[i] = timeStep*dy[i];
00199         myk6[i] = rCurrentYValues[i] - m8o27*mk1[i] + 2*mk2[i] - m3544o2565*mk3[i]
00200                         + m1859o4104*mk4[i] - 0.275*mk5[i];
00201     }
00202 
00203     pAbstractOdeSystem->EvaluateYDerivatives(time+0.5*timeStep, myk6, dy);
00204     for (unsigned i=0; i<num_equations; i++)
00205     {
00206         mk6[i] = timeStep*dy[i];
00207         mError[i] = (1/timeStep)*fabs(m1o360*mk1[i] - m128o4275*mk3[i]
00208                     - m2197o75240*mk4[i] + 0.02*mk5[i]+ m2o55*mk6[i]);
00209         rNextYValues[i] = rCurrentYValues[i] + m25o216*mk1[i] + m1408o2565*mk3[i]
00210                         + m2197o4104*mk4[i] - 0.2*mk5[i];
00211     }
00212 }
00213 
00214 void RungeKuttaFehlbergIvpOdeSolver::AdjustStepSize(double& rCurrentStepSize,
00215                                 const double& rError,
00216                                 const double& rTolerance,
00217                                 const double& rMaxTimeStep,
00218                                 const double& rMinTimeStep)
00219 {
00220     // Work out scaling factor delta for the step size
00221     double delta = pow(rTolerance/(2.0*rError), 0.25);
00222 
00223     // Maximum adjustment is *0.1 or *4
00224     if (delta <= 0.1)
00225     {
00226         rCurrentStepSize *= 0.1;
00227     }
00228     else if (delta >= 4.0)
00229     {
00230         rCurrentStepSize *= 4.0;
00231     }
00232     else
00233     {
00234         rCurrentStepSize *= delta;
00235     }
00236 
00237     if (rCurrentStepSize > rMaxTimeStep)
00238     {
00239         rCurrentStepSize = rMaxTimeStep;
00240     }
00241 
00242     if (rCurrentStepSize < rMinTimeStep)
00243     {
00244         std::cout << "rCurrentStepSize = " << rCurrentStepSize << "\n" << std::flush;
00245         std::cout << "rMinTimeStep = " << rMinTimeStep << "\n" << std::flush;
00246 
00247         EXCEPTION("RKF45 Solver: Ode needs a smaller timestep than the set minimum\n");
00248     }
00249 
00250 }
00251 
00252 
00253 /*
00254  * PUBLIC FUNCTIONS=============================================================
00255  */
00256 
00257 RungeKuttaFehlbergIvpOdeSolver::RungeKuttaFehlbergIvpOdeSolver()
00258     : m1932o2197(1932.0/2197.0), // you need the .0 s - caused me no end of trouble.
00259       m7200o2197(7200.0/2197.0),
00260       m7296o2197(7296.0/2197.0),
00261       m12o13(12.0/13.0),
00262       m439o216(439.0/216.0),
00263       m3680o513(3680.0/513.0),
00264       m845o4104(845.0/4104.0),
00265       m8o27(8.0/27.0),
00266       m3544o2565(3544.0/2565.0),
00267       m1859o4104(1859.0/4104.0),
00268       m1o360(1.0/360.0),
00269       m128o4275(128.0/4275.0),
00270       m2197o75240(2197.0/75240.0),
00271       m2o55(2.0/55.0),
00272       m25o216(25.0/216.0),
00273       m1408o2565(1408.0/2565.0),
00274       m2197o4104(2197.0/4104.0)
00275 {
00276 }
00277 
00278 OdeSolution RungeKuttaFehlbergIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00279                                                   std::vector<double>& rYValues,
00280                                                   double startTime,
00281                                                   double endTime,
00282                                                   double timeStep,
00283                                                   double tolerance)
00284 {
00285     assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00286     assert(endTime > startTime);
00287     assert(timeStep > 0.0);
00288 
00289     mStoppingEventOccurred = false;
00290     if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00291     {
00292         EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
00293     }
00294     // Allocate working memory
00295     std::vector<double> working_memory(rYValues.size());
00296     // And solve...
00297     OdeSolution solutions;
00298     //solutions.SetNumberOfTimeSteps((unsigned)(10.0*(startTime-endTime)/timeStep));
00299     bool return_solution = true;
00300     InternalSolve(solutions, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-5, tolerance, return_solution);
00301     return solutions;
00302 }
00303 
00304 void RungeKuttaFehlbergIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00305                                            std::vector<double>& rYValues,
00306                                            double startTime,
00307                                            double endTime,
00308                                            double timeStep)
00309 {
00310     assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00311     assert(endTime > startTime);
00312     assert(timeStep > 0.0);
00313 
00314     mStoppingEventOccurred = false;
00315     if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00316     {
00317         EXCEPTION("(Solve without sampling) Stopping event is true for initial condition");
00318     }
00319     // Allocate working memory
00320     std::vector<double> working_memory(rYValues.size());
00321     // And solve...
00322     OdeSolution not_required_solution;
00323     bool return_solution = false;
00324     InternalSolve(not_required_solution, pOdeSystem, rYValues, working_memory, startTime, endTime, timeStep, 1e-4, 1e-5, return_solution);
00325 }

Generated on Tue Aug 4 16:10:24 2009 for Chaste by  doxygen 1.5.5