Page MenuHomec4science
Files F86105787

AdamsMoulton.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Oct 4, 07:43

AdamsMoulton.cpp
View Options

#include <Eigen/Dense>
#include "AdamsMoulton.h"
AdamsMoulton::AdamsMoulton():MultistepAbstractSolver(){
};
AdamsMoulton::AdamsMoulton(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
convergenceRate=steps+1;
};
void AdamsMoulton::setStepsAlgo(int n){
stepsAlgo=n;
convergenceRate=n+1;
setBCoeff();
}
void AdamsMoulton::setBCoeff() {
bCoefficients.resize(stepsAlgo+1);
switch(stepsAlgo){
case 1:
bCoefficients[0]=0.5,bCoefficients[1]=0.5;
break;
case 2:
bCoefficients[0]= -1./12.,bCoefficients[1]=8./12.,bCoefficients[2]=5./12.;
break;
case 3:
bCoefficients[0]=1./24.,bCoefficients[1]=-5./24.,bCoefficients[2]=19./24.,bCoefficients[3]= 9./24.;
break;
case 4:
bCoefficients[0]=-19./720.,bCoefficients[1]=106./720.,bCoefficients[2]= -264./720.,bCoefficients[3]=646./720.,bCoefficients[4]=251./720.;
break;
default:
throw std::runtime_error("The number of steps of Adams Moulton has to be an integer between 1 and 4.");
}
}
void AdamsMoulton::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
// Method of AbstractSolver mother class is called.
AbstractSolver::setRightHandSide(rhsMatrix, *g_);
// Matrix I is needed to compute matrix defining iterations of Adams Moulton.
Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
// The coefficients defining the method are set.
setBCoeff();
//The matrix updateMatrix involved in the computation of the next step value of y is computed.
int s=bCoefficients.size();
double d=bCoefficients[s-1]*stepSize;
updateMatrix=I-d*A;
};
void AdamsMoulton::setRightHandSide(Eigen::MatrixXd & rhsMatrix, std::function<Eigen::VectorXd(double)> g_ ){
AbstractSolver::setRightHandSide(rhsMatrix, g_);
Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
setBCoeff();
int s=bCoefficients.size();
double d=bCoefficients[s-1]*stepSize;
updateMatrix=I-d*A;
};
void AdamsMoulton::step(Eigen::VectorXd &y, double t){
// rhsVector of the system to be solved at the current step is computed
Eigen::VectorXd rhsVector(dimension);
int s=bCoefficients.size();
rhsVector=y+stepSize*bCoefficients[s-1]*g(t+stepSize);
int count=0;
auto it=rhsPreviousSteps.begin();
for(it=rhsPreviousSteps.begin();it!=rhsPreviousSteps.end();it++){
rhsVector+=(*it)*bCoefficients[count]*stepSize;
count++;
}
// Solution at next time step is computed.
Eigen::VectorXd yNext=updateMatrix.lu().solve(rhsVector);
// The value of the rhs function at the next time step is computed and stored in rhsPreviousSteps.
Eigen::VectorXd rhsNext=A*yNext+g(t+stepSize);
rhsPreviousSteps.push_back(rhsNext);
// The first value of the list rhsPreviousSteps is not needed for the next integration step, so it is discarded.
rhsPreviousSteps.pop_front();
//update of y
y=yNext;
};

Event Timeline

c4science · Help