37#include "BackwardEulerIvpOdeSolver.hpp"
43 std::vector<double>& rCurrentYValues,
44 std::vector<double>& rCurrentGuess)
50 mResidual[i] = rCurrentGuess[i] - timeStep * dy[i] - rCurrentYValues[i];
57 std::vector<double>& rCurrentYValues,
58 std::vector<double>& rCurrentGuess)
117 if (fabs(pVector[i]) > norm)
119 norm = fabs(pVector[i]);
128 std::vector<double>& rCurrentYValues,
129 std::vector<double>& rCurrentGuess)
137 ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, rCurrentGuess);
143 for (
unsigned global_column=0; global_column<
mSizeOfOdeSystem; global_column++)
147 guess_perturbed[i] = rCurrentGuess[i];
150 guess_perturbed[global_column] += epsilon;
152 ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, guess_perturbed);
159 double one_over_eps = 1.0/epsilon;
162 mJacobian[i][global_column] = one_over_eps*(residual_perturbed[i] - residual[i]);
170 std::vector<double>& rCurrentYValues,
171 std::vector<double>& rNextYValues)
176 const double eps = 1e-6;
180 current_guess.assign(rCurrentYValues.begin(), rCurrentYValues.end());
184 unsigned counter = 0u;
185 const unsigned iter_limit = 20u;
186 while (norm > eps && counter < iter_limit)
189 ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
190 ComputeJacobian(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
201 current_guess[i] -=
mUpdate[i];
206 assert(counter < iter_limit);
208 rNextYValues.assign(current_guess.begin(), current_guess.end());
#define CHASTE_CLASS_EXPORT(T)
virtual void AnalyticJacobian(const std::vector< double > &rSolutionGuess, double **jacobian, double time, double timeStep)=0
virtual void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)=0
bool GetUseAnalyticJacobian()
unsigned GetNumberOfStateVariables() const
void ForceUseOfNumericalJacobian()
void ComputeNumericalJacobian(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rCurrentGuess)
double mNumericalJacobianEpsilon
void CalculateNextYValue(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rNextYValues)
unsigned mSizeOfOdeSystem
~BackwardEulerIvpOdeSolver()
void ComputeResidual(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rCurrentGuess)
void SetEpsilonForNumericalJacobian(double epsilon)
void ComputeJacobian(AbstractOdeSystem *pAbstractOdeSystem, double timeStep, double time, std::vector< double > &rCurrentYValues, std::vector< double > &rCurrentGuess)
bool mForceUseOfNumericalJacobian
BackwardEulerIvpOdeSolver(unsigned sizeOfOdeSystem)
double ComputeNorm(double *pVector)