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
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);
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
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
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
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
00167 assert(mSizeOfOdeSystem == pAbstractOdeSystem->GetNumberOfStateVariables());
00168
00169 unsigned counter = 0;
00170
00171 const double eps = 1e-6;
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
00180 ComputeResidual(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
00181 ComputeJacobian(pAbstractOdeSystem, timeStep, time, rCurrentYValues, current_guess);
00182
00183
00184
00185
00186 SolveLinearSystem();
00187
00188
00189 norm = ComputeNorm(mUpdate);
00190
00191
00192 for (unsigned i=0; i<mSizeOfOdeSystem; i++)
00193 {
00194 current_guess[i] -= mUpdate[i];
00195 }
00196
00197 counter++;
00198 assert(counter < 20);
00199 }
00200 rNextYValues.assign(current_guess.begin(), current_guess.end());
00201 }
00202
00203 BackwardEulerIvpOdeSolver::BackwardEulerIvpOdeSolver(unsigned sizeOfOdeSystem)
00204 {
00205 mSizeOfOdeSystem = sizeOfOdeSystem;
00206
00207
00208 mNumericalJacobianEpsilon = 1e-6;
00209 mForceUseOfNumericalJacobian = false;
00210
00211
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
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
00248 #include "SerializationExportWrapperForCpp.hpp"
00249 CHASTE_CLASS_EXPORT(BackwardEulerIvpOdeSolver)