Page MenuHomec4science
Files F63618266

main.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, May 21, 09:43

main.cpp
View Options

#include <iostream>
#include <fstream>
#include <Eigen/Dense>
#include <vector>
#include <cmath>
#include "ODElibrary/ForwardEuler.h"
#include "ODElibrary/RungeKutta.h"
#include "ODElibrary/AdamsBashforth.h"
#include "ODElibrary/AdamsMoulton.h"
#include "ODElibrary/BDF.h"
/**
* This function will correspond to the rhs function g(t), where:
* @param t: time at which it has to be computed.
* It is designed to compute the values at a single specific time.
*/
Eigen::VectorXd g_function(double t);
/**\brief This function corresponds to another possible rhs function g(t) (more details below).
*/
Eigen::VectorXd ForcedOscillations(double t);
int main(int argc, char *argv[]){
std::cout << "Number of command line arguments = "
<< argc << "\n";
for (int i=0; i<argc; i++) {
std::cout << "Argument " << i << " = " << argv[i];
std::cout << "\n";
}
Eigen::VectorXd initialValues;
double parameters[11];
std::fstream inputfile;
inputfile.open(argv[1]);
std:: string s;
for (int i=0; i<11;i++) {
inputfile >> s >> parameters[i];
if (i==7){
int dimension=(int)(parameters[6]);
initialValues.resize(dimension);
initialValues[0]=parameters[7];
int j=1;
while(j<dimension){
inputfile>>initialValues[j];
j++;
}
}
}
inputfile.close();
std::vector<AbstractSolver*> solvers;
int solver_name=(int)(parameters[0]);
switch(solver_name){
case 1: { //Euler
auto solver = new ForwardEuler();
solvers.push_back(solver);
break;
}
case 2: { //RK
auto solver = new RungeKutta();
solver->setMethodOrder((int)(parameters[10]));
solvers.push_back(solver);
break;
}
case 3: { //Adams Bashforth
auto solver = new AdamsBashforth();
solver->setStepsAlgo((int)(parameters[9]));
solvers.push_back(solver);
break;
}
case 4: { // Adams Moulton
auto solver = new AdamsMoulton();
solver->setStepsAlgo((int)(parameters[9]));
solvers.push_back(solver);
break;
}
case 5: { //BDF
auto solver = new BDF();
solver->setStepsAlgo((int)(parameters[9]));
solvers.push_back(solver);
break;
}
default:
std::cerr<<"Method not known. Set a number from 1 to 5.\n";
}
AbstractSolver* solver=solvers[0];
solver->setInitialTime(parameters[1]);
solver->setFinalTime(parameters[2]);
int discretizationMethod=parameters[5];
if(discretizationMethod==0)
solver->setNumberOfSteps(parameters[3]);
else
solver->setStepSize(parameters[4]);
solver->setInitialValue(initialValues);
Eigen::VectorXd (*p_function)(double t);
Eigen::MatrixXd m(2,2);
if((int)(parameters[8])==1) {
p_function = &g_function;
m(0,0)=1.0, m(0,1)=0.0,m(1,0)=0.0, m(1,1)=-2.0;
}
else{
p_function=&ForcedOscillations;
double springStiffness=-0.001;
m(0,0)=0, m(0,1)=1,m(1,0)=-springStiffness,m(1,1)=0.3;
}
solver->setRightHandSide(m,p_function);
std::ofstream outputfile;
outputfile.open("solution.dat");
solver->solve(outputfile);
outputfile.close();
/*
Eigen::VectorXd y0(2);
y0[0]=1.,y0[1]=1.;
/**
* Another possible problem that can be used to test this library describes the behaviour of a spring with a mass attached to it,
* with an additional term representing forced oscillations. The matrix describing it is here defined, along with the
* initial vector used to solve it, y_1(0)=0, y_2(0)=0. At t=1, the exact solutions of the vectorial differential equation
* is y_1(1)=0.182532458, y_2(1)=0.046017418.
double springStiffness=0.1;
Eigen::MatrixXd p(2,2); // from spring equation
p(0,0)=0, p(0,1)=1,p(1,0)=-springStiffness,p(1,1)=0.3;
Eigen::VectorXd (*p_forcedoscillations)(double t);
p_forcedoscillations=&ForcedOscillations;
Eigen::VectorXd yStart(2);
yStart[0]=0.,yStart[1]=0.;
/**
* In the following part of the code, different solvers can be called to solve the ODE problem at hand.
* By leaving uncommented -ONLY- the line corresponding to the method which one wants to apply, the problem
* is solved with it. Parameters can be played with.
double initialTime=0;
double finalTime=1;
int nSteps=100;
int nStepsMultistep=5;
int orderRK=2;
/**\brief Change y0 to yStart to integrate the spring problem instead of the exponential test problem.
//ForwardEuler solverq=ForwardEuler(initialTime,finalTime,nSteps,yStart);
//RungeKutta solverq=RungeKutta(initialTime,finalTime,nSteps,orderRK,yStart);
//AdamsBashforth solverq=AdamsBashforth(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
//AdamsMoulton solverq=AdamsMoulton(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
//BDF solverq=BDF(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
/**
* \brief The rhs function is set to the desired value. Any matrix and any pointer to a function returning a
* Eigen::VectorXd can be passed as parameters. The solver will then use them to integrate the differential equations.
* Setting solver.setRightHandSide(p,p_forcedoscillations), the spring problem described above can be solved.
//solverq.setRightHandSide(m,p_function);
//solverq.setRightHandSide(p,p_forcedoscillations);
/**
* The problem is solved and the solution is displayed on std::cout, since this is the default of the 'solve' method.
* By calling solver.solve(solutionfile) the computed values will instead be written to the file 'solutionfile'.
std::ofstream solutionfile;
solutionfile.open("AdamsBashforth2_solution.txt");
//solverq.solve();
solutionfile.close();
*/
}
/**
* Implementation of the rhs functions g. The entries of the result are set to zero so that the accuracy of the simple integration
* performed in this main can be easily checked. It can be changed to whatever is desired by the user. The return value
* does not have to be two-dimensional, the function can be modified to return a vector of whatever length.
* It is set to two dimensions just for demonstration purposes.
* @param t : time at which the function is computed
* @return : An Eigen::VectorXd storing the value of g(t).
*/
Eigen::VectorXd g_function(double t){
Eigen::VectorXd result(2);
result[0]=0;
result[1]=0;
return result;
}
Eigen::VectorXd ForcedOscillations(double t){
Eigen::VectorXd result(2);
result[0] = 0;
result[1] = cos(2*M_PI*t); // derivative of velocity
return result;
};

Event Timeline

c4science · Help