BackwardEulerIvpOdeSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 
00030 #include "BackwardEulerIvpOdeSolver.hpp"
00031 #include <cmath>
00032 
00033 void BackwardEulerIvpOdeSolver::ComputeResidual(AbstractOdeSystem* pAbstractOdeSystem,
00034                                                 double timeStep,
00035                                                 double time,
00036                                                 std::vector<double>& rCurrentYValues,
00037                                                 std::vector<double>& rCurrentGuess)
00038 {
00039     std::vector<double> dy(mSizeOfOdeSystem); //For JC to optimize
00040     pAbstractOdeSystem->EvaluateYDerivatives(time+timeStep, rCurrentGuess, dy);
00041     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00042     {
00043         mResidual[i] = rCurrentGuess[i] - timeStep * dy[i] - rCurrentYValues[i];
00044     }
00045 }
00046 
00047 void BackwardEulerIvpOdeSolver::ComputeJacobian(AbstractOdeSystem* pAbstractOdeSystem,
00048                                                 double timeStep,
00049                                                 double time,
00050                                                 std::vector<double>& rCurrentYValues,
00051                                                 std::vector<double>& rCurrentGuess)
00052 {
00053     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00054     {
00055         for (unsigned j=0; j<mSizeOfOdeSystem; j++)
00056         {
00057             mJacobian[i][j] = 0.0;
00058         }
00059     }
00060 
00061     if (pAbstractOdeSystem->GetUseAnalyticJacobian() && !mForceUseOfNumericalJacobian)
00062     {
00063         // The ODE system has an analytic jacobian, so use that
00064         AbstractOdeSystemWithAnalyticJacobian* p_ode_system
00065             = static_cast<AbstractOdeSystemWithAnalyticJacobian*>(pAbstractOdeSystem);
00066         p_ode_system->AnalyticJacobian(rCurrentGuess, mJacobian, time, timeStep);
00067     }
00068     else
00069     {
00070         ComputeNumericalJacobian(pAbstractOdeSystem,
00071                                  timeStep,
00072                                  time,
00073                                  rCurrentYValues,
00074                                  rCurrentGuess);
00075     }
00076 }
00077 
00078 void BackwardEulerIvpOdeSolver::SolveLinearSystem()
00079 {
00080     double fact;
00081     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00082     {
00083         for (unsigned ii=i+1; ii<mSizeOfOdeSystem; ii++)
00084         {
00085             fact = mJacobian[ii][i]/mJacobian[i][i];
00086             for (unsigned j=i; j<mSizeOfOdeSystem; j++)
00087             {
00088                 mJacobian[ii][j] -= fact*mJacobian[i][j];
00089             }
00090             mResidual[ii] -= fact*mResidual[i];
00091         }
00092     }
00093     // This needs to int, since a downloop in unsigned won't terminate properly
00094     for (int i=mSizeOfOdeSystem-1; i>=0; i--)
00095     {
00096         mUpdate[i] = mResidual[i];
00097         for (unsigned j=i+1; j<mSizeOfOdeSystem; j++)
00098         {
00099             mUpdate[i] -= mJacobian[i][j]*mUpdate[j];
00100         }
00101         mUpdate[i] /= mJacobian[i][i];
00102     }
00103 }
00104 
00105 double BackwardEulerIvpOdeSolver::ComputeNorm(double* pVector)
00106 {
00107     double norm = 0.0;
00108     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00109     {
00110         if (fabs(pVector[i]) > norm)
00111         {
00112             norm = fabs(pVector[i]);
00113         }
00114     }
00115     return norm;
00116 }
00117 
00118 void BackwardEulerIvpOdeSolver::ComputeNumericalJacobian(AbstractOdeSystem* pAbstractOdeSystem,
00119                                                          double timeStep,
00120                                                          double time,
00121                                                          std::vector<double>& rCurrentYValues,
00122                                                          std::vector<double>& rCurrentGuess)
00123 {
00124     std::vector<double> residual(mSizeOfOdeSystem);
00125     std::vector<double> residual_perturbed(mSizeOfOdeSystem);
00126     std::vector<double> guess_perturbed(mSizeOfOdeSystem);
00127 
00128     double epsilon = mNumericalJacobianEpsilon;
00129 
00130     ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, rCurrentGuess);
00131     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00132     {
00133         residual[i] = mResidual[i];
00134     }
00135 
00136     for (unsigned global_column=0; global_column<mSizeOfOdeSystem; global_column++)
00137     {
00138         for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00139         {
00140             guess_perturbed[i] = rCurrentGuess[i];
00141         }
00142 
00143         guess_perturbed[global_column] += epsilon;
00144 
00145         ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, guess_perturbed);
00146         for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00147         {
00148             residual_perturbed[i] = mResidual[i];
00149         }
00150 
00151         // Compute residual_perturbed - residual
00152         double one_over_eps = 1.0/epsilon;
00153         for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00154         {
00155             mJacobian[i][global_column] = one_over_eps*(residual_perturbed[i] - residual[i]);
00156         }
00157     }
00158 }
00159 
00160 void BackwardEulerIvpOdeSolver::CalculateNextYValue(AbstractOdeSystem* pAbstractOdeSystem,
00161                                                     double timeStep,
00162                                                     double time,
00163                                                     std::vector<double>& rCurrentYValues,
00164                                                     std::vector<double>& rNextYValues)
00165 {
00166     // check the size of the ode system matches the solvers expected
00167     assert(mSizeOfOdeSystem == pAbstractOdeSystem->GetNumberOfStateVariables());
00168 
00169     unsigned counter = 0;
00170 //        const double eps = 1e-6 * rCurrentGuess[0]; // Our tolerance (should use min(guess) perhaps?)
00171     const double eps = 1e-6; // JonW tolerance
00172     double norm = 2*eps;
00173 
00174     std::vector<double> current_guess(mSizeOfOdeSystem);
00175     current_guess.assign(rCurrentYValues.begin(), rCurrentYValues.end());
00176 
00177     while (norm > eps)
00178     {
00179         // Calculate Jacobian and mResidual for current guess
00180         ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
00181         ComputeJacobian(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
00182 //            // Update norm (our style)
00183 //            norm = ComputeNorm(mResidual);
00184 
00185         // Solve Newton linear system
00186         SolveLinearSystem();
00187 
00188         // Update norm (JonW style)
00189         norm = ComputeNorm(mUpdate);
00190 
00191         // Update current guess
00192         for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00193         {
00194             current_guess[i] -= mUpdate[i];
00195         }
00196 
00197         counter++;
00198         assert(counter < 20); // avoid infinite loops
00199     }
00200     rNextYValues.assign(current_guess.begin(), current_guess.end());
00201 }
00202 
00203 BackwardEulerIvpOdeSolver::BackwardEulerIvpOdeSolver(unsigned sizeOfOdeSystem)
00204 {
00205     mSizeOfOdeSystem = sizeOfOdeSystem;
00206 
00207     // default epsilon
00208     mNumericalJacobianEpsilon = 1e-6;
00209     mForceUseOfNumericalJacobian = false;
00210 
00211     // allocate memory
00212     mResidual = new double[mSizeOfOdeSystem];
00213     mUpdate = new double[mSizeOfOdeSystem];
00214 
00215     mJacobian = new double*[mSizeOfOdeSystem];
00216     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00217     {
00218         mJacobian[i] = new double[mSizeOfOdeSystem];
00219     }
00220 }
00221 
00222 BackwardEulerIvpOdeSolver::~BackwardEulerIvpOdeSolver()
00223 {
00224     // Delete pointers
00225     delete[] mResidual;
00226     delete[] mUpdate;
00227 
00228     for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00229     {
00230         delete[] mJacobian[i];
00231     }
00232     delete[] mJacobian;
00233 }
00234 
00235 void BackwardEulerIvpOdeSolver::SetEpsilonForNumericalJacobian(double epsilon)
00236 {
00237     assert(epsilon > 0);
00238     mNumericalJacobianEpsilon = epsilon;
00239 }
00240 
00241 void BackwardEulerIvpOdeSolver::ForceUseOfNumericalJacobian()
00242 {
00243     mForceUseOfNumericalJacobian = true;
00244 }
00245 
00246 
00247 // Serialization for Boost >= 1.36
00248 #include "SerializationExportWrapperForCpp.hpp"
00249 CHASTE_CLASS_EXPORT(BackwardEulerIvpOdeSolver)

Generated on Mon Nov 1 12:35:24 2010 for Chaste by  doxygen 1.5.5