Homec4science
Diffusion PCSCproject (master) Browse

R9686/ODElibrary/BDF.cpp
5a4a7bdba305master

Last Change
Raw File

BDF.cpp
View Options

#include <Eigen/Dense>
#include "BDF.h"
BDF::BDF():MultistepAbstractSolver(){
};
BDF::BDF(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
convergenceRate=steps;
};
void BDF::setStepsAlgo(int n){
stepsAlgo=n;
convergenceRate=n;
setBCoeff();
}
void BDF::setBCoeff() {
bCoefficients.resize(stepsAlgo+1);
switch(stepsAlgo){
case 1:
bCoefficients[0]=-1.,bCoefficients[1]=1.;
break;
case 2:
bCoefficients[0]= 0.5,bCoefficients[1]=-2.,bCoefficients[2]=1.5;
break;
case 3:
bCoefficients[0]=-1./3.,bCoefficients[1]=1.5,bCoefficients[2]= -3.,bCoefficients[3]=11./6.;
break;
case 4:
bCoefficients[0]=1./4.,bCoefficients[1]=-4./3.,bCoefficients[2]= 3.,bCoefficients[3]=-4.,bCoefficients[4]= 25./12.;
break;
case 5:
bCoefficients[0]=-1./5.,bCoefficients[1]=5./4.,bCoefficients[2]= -10./3.,bCoefficients[3]=5.,bCoefficients[4]=-5.,bCoefficients[5]=137./60.;
break;
default:
throw std::runtime_error("The number of steps of Adams Moulton has to be an integer between 1 and 5.");
}
}
void BDF::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 BDF.
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];
updateMatrix=d*I-stepSize*A;
};
void BDF::setRightHandSide(Eigen::MatrixXd & rhsMatrix, std::function<Eigen::VectorXd(double)> g_ ){
// Method of AbstractSolver mother class is called.
AbstractSolver::setRightHandSide(rhsMatrix, g_);
// Matrix I is needed to compute matrix defining iterations of BDF.
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];
updateMatrix=d*I-stepSize*A;
};
void BDF::step(Eigen::VectorXd &y, double t){
// rhsVector of the system to be solved at the current step is computed
Eigen::VectorXd rhsVector(dimension);
rhsVector=stepSize*g(t+stepSize);
int count=0;
//The vector yInitialConditions originally contains the missing initial conditions computed in order to initialize the algorithm.
//In the BDF method, this vector is then reused to store the last values of y, needed to perform the next step update.
auto it=yInitialConditions.begin();
for(it=yInitialConditions.begin();it!=yInitialConditions.end();it++){
rhsVector-=(*it)*bCoefficients[count];
count++;
}
// Solution at next time step is computed.
Eigen::VectorXd yNext=updateMatrix.lu().solve(rhsVector);
y=yNext;
yInitialConditions.push_back(yNext);
//The first value of the list is not needed for the next integration step, so it is discarded.
yInitialConditions.pop_front();
};
c4science · Help