ode: tutorial 1#
#include "pressio/ode_advancers.hpp"
#include "pressio/ode_steppers_explicit.hpp"
#include <Eigen/Core>
template<class ScalarType>
class SimpleSystem{
public:
using independent_variable_type = ScalarType;
using state_type = Eigen::Matrix<ScalarType,-1,1>;
using rhs_type = state_type;
state_type createState() const{
auto s = state_type(3);
s.setConstant(0);
return s;
};
rhs_type createRhs() const{
auto v = rhs_type(3);
v.setConstant(0);
return v;
};
void rhs(const state_type & state,
const independent_variable_type timeIn,
rhs_type & rhs) const
{
constexpr ScalarType ten{10};
rhs(0) = ten * state(0);
rhs(1) = ten * state(1);
rhs(2) = ten * state(2);
};
};
int main(int argc, char *argv[])
{
pressio::log::initialize(pressio::logto::terminal);
using scalar_type = double;
SimpleSystem<scalar_type> problem;
auto y = problem.createState();
namespace pode = pressio::ode;
constexpr auto scheme = pode::StepScheme::ForwardEuler;
auto stepperObj = pode::create_explicit_stepper(scheme, problem);
{
/*
run a fixed number of steps
*/
y(0) = 1.; y(1) = 2.; y(2) = 3.;
const scalar_type dt{0.1};
const scalar_type startTime{0.0};
pode::advance_n_steps(stepperObj, y, startTime, dt, pode::StepCount{1});
// check solution
std::cout << "Computed solution: ["
<< y(0) << " " << y(1) << " " << y(2) << "] "
<< "Expected solution: [2,4,6] \n";
}
std::cout << std::endl;
{
/*
run until target time
*/
y(0) = 1.; y(1) = 2.; y(2) = 3.;
const scalar_type startTime{0.0};
const scalar_type finalTime{1.0};
pode::advance_to_target_point(stepperObj, y, startTime, finalTime,
[](const pode::StepCount & /*unused*/,
const pode::StepStartAt<scalar_type> & /*unused*/,
pode::StepSize<scalar_type> & dt)
{
dt = 0.1;
});
// check solution
std::cout << "Computed solution: ["
<< y(0) << " " << y(1) << " " << y(2) << "] "
<< "Expected solution: [1024, 2048, 3074] \n";
}
pressio::log::finalize();
return 0;
}