Page MenuHomec4science
Files F63472874

MultistepAbstractSolver.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 20, 09:17

MultistepAbstractSolver.cpp
View Options

#include "MultistepAbstractSolver.h"
MultistepAbstractSolver::MultistepAbstractSolver(){
initialTime = 0;
finalTime = 0;
nSteps = 0;
stepSize = 0;
stepsAlgo=0;
dimension=0;
convergenceRate=0;
};
MultistepAbstractSolver::MultistepAbstractSolver(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):AbstractSolver(t0,tf,n,y0) {
stepsAlgo=steps;
yInitialConditions.push_back(y0);
};
int MultistepAbstractSolver::getStepsAlgo(){
return stepsAlgo;
};
void MultistepAbstractSolver::computeMissingInitialValues(int orderRK){
//The coefficients for computing missing initial conditions through Runge Kutta are initialized.
Eigen::Vector4d aCoeff;
switch(orderRK){
case 0:
{
rhsPreviousSteps.clear();
Eigen::VectorXd rhsStart=A*initialValue+g(initialTime);
rhsPreviousSteps.push_back(rhsStart);
return;
}
break;
case 1:
aCoeff[0]=1, aCoeff[1]=0, aCoeff[2]=0,aCoeff[3]=0;
break;
case 2:
aCoeff[0]=0.5, aCoeff[1]=0.5, aCoeff[2]=0,aCoeff[3]=0;
break;
case 3:
aCoeff[0]=1./6., aCoeff[1]=4./6., aCoeff[2]=1./6.,aCoeff[3]=0;
break;
case 4:
aCoeff[0]=1./4., aCoeff[1]=1./4., aCoeff[2]=1./4.,aCoeff[3]=1./4.;
break;
default:
throw std::runtime_error("Order orderRK required too high. Try using a multistep method with fewer steps.");
}
auto yPrev=initialValue;
auto tPrev=initialTime;
//rhsPreviousSteps list is cleared and first value is stored. Need of clearing it comes from test function.
rhsPreviousSteps.clear();
Eigen::VectorXd rhsStart=A*initialValue+g(initialTime);
rhsPreviousSteps.push_back(rhsStart);
// k1, k2, k3, k4 define the updates of RK
Eigen::VectorXd k1(dimension),k2(dimension),k3(dimension),k4(dimension);
k2.setZero(), k3.setZero(),k4.setZero();
int initialConditionsNeeded=orderRK;
if(stepsAlgo==convergenceRate-1)
initialConditionsNeeded-=1;
for (int i=0;i<initialConditionsNeeded;i++){
k1=A*yPrev+g(tPrev);
//k2, k3, k4 are computed only if a higher order RK method is needed.
if(orderRK>1)
k2=A*(yPrev+0.5*stepSize*k1)+g(tPrev+0.5*stepSize); // With this choice, the order 2 method corresponds to Euler's modified method.
if(orderRK>2)
k3 = A * (yPrev + stepSize * (-1 * k1 + 2 * k2)) +g(tPrev + stepSize);// This choice for k3 allows us to "recycle" k1, k2, obtaining a 3rd order RK.
if(orderRK>3) {
k3=A *(yPrev+0.5*stepSize*k2)+g(tPrev+0.5*stepSize);
k4 = A * (yPrev + stepSize*k3) + g(tPrev + stepSize);
}
//computation of a missing initial condition.
Eigen::VectorXd yNext=yPrev+stepSize*(aCoeff[0]*k1+aCoeff[1]*k2+aCoeff[2]*k3+aCoeff[3]*k4);
//rhs function evaluated and stored in rhsPreviousSteps.
Eigen::VectorXd rhsNext=A*yNext+g(tPrev+stepSize);
rhsPreviousSteps.push_back(rhsNext);
//The computed missing initial condition is stored in the list yInitialConditions.
yInitialConditions.push_back(yNext);
yPrev=yNext;
tPrev+=stepSize;
}
}
void MultistepAbstractSolver::solve(std::ostream &stream){
yInitialConditions.clear();
yInitialConditions.push_back(initialValue);
// The missing initial values are computed, with a suitable one step algorithm.
if(stepsAlgo==convergenceRate-1)
// For implicit method AdamsMoulton, a higher order RK method is needed to compute the missing initial values.
computeMissingInitialValues(convergenceRate-1);
else
// For AdamsBashforth, BDF, stepsAlgo==convergenceRate therefore a lower RK is needed w.r.t. AdamsMoulton.
computeMissingInitialValues(stepsAlgo-1);
//The initial conditions computed with computeMissingInitialValues are streamed to the output stream chosen by the user.
int count=0;
auto it=yInitialConditions.begin();
for(it=yInitialConditions.begin();it!=yInitialConditions.end();it++){
stream<<initialTime+stepSize*count<<" "<<(*it).transpose()<<"\n";
count++;
}
// t is the time at which y has to be computed at each step; y is initialized with the value of the last missing initial
// condition computed; its value is progressively changed by the algorithm's iterations.
double t=initialTime+stepSize*(stepsAlgo-1);
Eigen::VectorXd y = yInitialConditions.back();
// nSteps-stepsAlgo steps are needed because missing initial conditions have already been computed.
for (int i=0;i<=nSteps-stepsAlgo;i++) {
step(y,t);
stream<<initialTime+stepSize*(i+stepsAlgo)<<" "<<y.transpose()<<"\n";
t+=stepSize;
}
}
//void MultistepAbstractSolver::step(Eigen::VectorXd &y, double t) {
// throw std::runtime_error("function step not implemented in MultistepAbstractSolver: it's a purely virtual method, implemented in the daughter classes.");
//}

Event Timeline

c4science · Help