Я пытаюсь сравнить ODE-решатель GSL с odeint from boost для многомерных систем и поэтому написал короткую тестовую программу:
#include <iostream>
#include <cmath>
#include <vector>
#include <chrono>
#include <boost/numeric/odeint.hpp>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
constexpr int bench_rounds = 10000;
constexpr int dim = 100;
typedef std::vector<double> state_type;
struct simulation_parameters{
int dimension = dim; /* number of differential equations */
double eps_abs = 1.e-8; /* absolute error requested */
double eps_rel = 1.e-11; /* relative error requested */
double tmin = 0.; /* starting t value */
double tmax = 10.; /* final t value */
double delta_t = 1e-2;
double h = 1e-6; /* starting step size for ode solver */
};
struct push_back_state_and_time
{
std::vector< state_type >& m_states;
std::vector< double >& m_times;
push_back_state_and_time( std::vector< state_type > &states , std::vector< double > × )
: m_states( states ) , m_times( times ) { }
void operator()( const state_type &x , double t )
{
m_states.push_back( x );
m_times.push_back( t );
}
};
int gsl_rhs (double t, const double y[], double f[], void *params_ptr)
{
/* get parameter(s) from params_ptr; here, just a double */
double mu = *(double *) params_ptr;
/* evaluate the right-hand-side functions at t */
for(size_t i = 0; i < dim; ++i)
f[i] = -2 * M_PI * y[i];
return GSL_SUCCESS; /* GSL_SUCCESS defined in gsl/errno.h as 0 */
};
void boost_rhs(const state_type &y, state_type &f, const int t){
(void) t;
for(size_t i = 0; i < y.size(); ++i)
f[i] = -2 * M_PI * y[i];
}
class boost_calculator{
public:
boost_calculator(const size_t n_dimensions) : n_dimensions(n_dimensions) {};
void operator()(const state_type &x, state_type &dx, const double){
dx.resize(n_dimensions);
for(size_t i = 0; i < n_dimensions; ++i)
dx[i] = -2 * M_PI * x[i];
}
private:
size_t n_dimensions;
};
double exact_solution(const double t){
return exp(-2 * M_PI * t);
};
double check_solution_correctness(const std::vector<std::pair<double, double>> &solution){
double error = 0;
for(size_t i = 0; i < solution.size(); ++i){
const double local_error = std::abs(solution[i].second - exact_solution (solution[i].first));
error += local_error;
}
return error;
}
double check_solution_correctness(const std::vector<std::pair<double, double*>> &solution){
double error = 0;
for(size_t i = 0; i < solution.size(); ++i){
double local_error = 0;
for(size_t j = 0; j < dim; ++j)
local_error += std::abs(solution[i].second[j] - exact_solution (solution[i].first));
error += local_error;
}
return error;
}
double check_solution_correctness(const std::vector<std::pair<double, std::vector<double>>> &solution){
double error = 0;
for(size_t i = 0; i < solution.size(); ++i){
double local_error = 0;
for(size_t j = 0; j < solution[i].second.size(); ++j)
local_error += std::abs(solution[i].second[j] - exact_solution (solution[i].first));
error += local_error;
}
return error;
}
double check_solution_correctness(const std::vector<std::pair<double, std::complex<double>>> &solution){
double error = 0;
for(size_t i = 0; i < solution.size(); ++i){
const double local_error = std::abs(solution[i].second - exact_solution (solution[i].first));
error += local_error;
}
return error;
}
void solve_equation_with_GSL(std::vector<std::pair<double, state_type>> &solution, const simulation_parameters ¶ms){
/* define the type of routine for making steps: */
const gsl_odeiv2_step_type *type_ptr = gsl_odeiv2_step_rkf45;
/*
allocate/initialize the stepper, the control function, and the
evolution function.
*/
gsl_odeiv2_step *step_ptr = gsl_odeiv2_step_alloc (type_ptr, params.dimension);
gsl_odeiv2_control *control_ptr = gsl_odeiv2_control_y_new (params.eps_abs, params.eps_rel);
gsl_odeiv2_evolve *evolve_ptr = gsl_odeiv2_evolve_alloc (params.dimension);
gsl_odeiv2_system my_system; /* structure with the rhs function, etc. */
double mu = 10; /* parameter for the diffeq */
double y[dim]; /* current solution vector */
double t, t_next; /* current and next independent variable */
double tmin, tmax, delta_t; /* range of t and step size for output */
double h = params.h;
tmin = params.tmin;
tmax = params.tmax;
delta_t = params.delta_t;
/* load values into the my_system structure */
my_system.function = gsl_rhs; /* the right-hand-side functions dy[i]/dt */
my_system.jacobian = nullptr; /* the Jacobian df[i]/dy[j] */
my_system.dimension = params.dimension; /* number of diffeq's */
my_system.params = μ /* parameters to pass to rhs and jacobian */
for(size_t i = 0; i < dim; ++i)
y[i] = 1.; /* initial x value */
t = tmin; /* initialize t */
//printf ("%.5e %.5e %.5e\n", t, y[0], exact_solution (t)); /* initial values */
/* step to tmax from tmin */
const int vector_size = (int)((tmax - tmin) / delta_t) + 1;
solution.clear();
solution.resize(vector_size - 1);
int counter = 0;
for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)
{
while (t < t_next) /* evolve from t to t_next */
{
gsl_odeiv2_evolve_apply (evolve_ptr, control_ptr, step_ptr,
&my_system, &t, t_next, &h, y);
}
//printf ("%.5e %.5e %.5e %.5e\n", t, y[0], y[1], exact_solution (t)); /* print at t=t_next */
solution[counter] = std::make_pair(t_next, state_type(y, y + dim));
counter++;
}
/* all done; free up the gsl_odeiv stuff */
gsl_odeiv2_evolve_free (evolve_ptr);
gsl_odeiv2_control_free (control_ptr);
gsl_odeiv2_step_free (step_ptr);
}
void solve_equation_with_Boost(std::vector<std::pair<double, state_type>> &solution, const simulation_parameters ¶ms){
using namespace boost::numeric::odeint;
typedef runge_kutta_fehlberg78<state_type> error_stepper_type;
state_type x_0 = state_type{1.};
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
std::vector<state_type> x_vec;
std::vector<double> t;
const size_t steps = integrate_adaptive(make_controlled<error_stepper_type>(params.eps_abs, params.eps_rel),
boost_rhs, x_0, params.tmin, params.tmax, params.delta_t, push_back_state_and_time(x_vec, t));
solution.clear();
solution.resize(steps);
for(size_t i = 0; i < steps; ++i)
solution[i] = std::make_pair(t[i], x_vec[i]);
}
void solve_equation_with_Boost_class(std::vector<std::pair<double, state_type>> &solution, const simulation_parameters ¶ms){
using namespace boost::numeric::odeint;
typedef runge_kutta_fehlberg78<state_type> error_stepper_type;
state_type x_0 = state_type{1.};
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
std::vector<state_type> x_vec;
std::vector<double> t;
boost_calculator local_calculator(dim);
const size_t steps = integrate_adaptive(make_controlled<error_stepper_type>(params.eps_abs, params.eps_rel),
local_calculator, x_0, params.tmin, params.tmax, params.delta_t, push_back_state_and_time(x_vec, t));
solution.clear();
solution.resize(steps);
for(size_t i = 0; i < steps; ++i)
solution[i] = std::make_pair(t[i], x_vec[i]);
}
int main()
{
std::vector<std::pair<double, state_type>> GSL_solution;
std::vector<std::pair<double, state_type>> boost_solution, boost_class_solution;
simulation_parameters local_sim_params;
auto t1 = std::chrono::high_resolution_clock::now();
for(size_t i = 0; i < bench_rounds; ++i)
solve_equation_with_GSL(GSL_solution, local_sim_params);
auto t2 = std::chrono::high_resolution_clock::now();
auto t3 = std::chrono::high_resolution_clock::now();
for(size_t i = 0; i < bench_rounds; ++i)
solve_equation_with_Boost (boost_solution, local_sim_params);
auto t4 = std::chrono::high_resolution_clock::now();
auto t5 = std::chrono::high_resolution_clock::now();
for(size_t i = 0; i < bench_rounds; ++i)
solve_equation_with_Boost_class (boost_class_solution, local_sim_params);
auto t6 = std::chrono::high_resolution_clock::now();
std::cout << "GSL-error: " << check_solution_correctness(GSL_solution)
<< ", solved within " << std::chrono::duration_cast<std::chrono::milliseconds>( t2 - t1 ).count()
<< " ms" << '\n';
std::cout << "Boost-error: " << check_solution_correctness (boost_solution)
<< ", solved within " << std::chrono::duration_cast<std::chrono::milliseconds>( t4 - t3 ).count()
<< " ms" << '\n';
std::cout << "Boost_class-error: " << check_solution_correctness (boost_class_solution)
<< ", solved within " << std::chrono::duration_cast<std::chrono::milliseconds>( t6 - t5 ).count()
<< " ms" << '\n';
return 0;
}
Изначально я начал использовать одно измерение и расширил это в нескольких измерениях. Мои текущие результаты теста:
GSL-error: 1.40024e-06, solved within 30924 ms
Boost-error: 1.74368e-08, solved within 1879 ms
Boost_class-error: 1.74368e-08, solved within 10482 ms
Насколько я знаю, мне нужно создать отдельный класс, если я хочу передать параметры в мою процедуру вычисления rhs, но сравнить ее с версией без классов, Я реализовал обе версии.
Теперь для dim = 1 этот код работал нормально, но для dim > 1 у меня возникли сбои в функции вычисления rhs-класса класса boost. Я мог исправить это, изменяя размер вектора решения каждый раз, когда я вызываю функцию, но это значительно замедляет эту функцию по сравнению с функцией без классов (как видно из результатов). Таким образом, есть ли другой способ реализации функции на основе классов без изменения размера вектора? И почему функции без классов всегда имеют правильный размер для f, а функции класса - нет?