将 boost::odeint 与向量类一起使用,而无需调整向量的大小
Use boost::odeint with class for vectors without resizing vectors
我正在尝试将GSL的ODE求解器与多维系统的boost的odeint进行比较,因此编写了一个简短的测试程序:
#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 %.5en", 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 %.5en", 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
,我在 boost 类的 rhs 计算函数中遇到了段错误。我可以修复它的方法是每次调用函数时调整解决方案向量的大小,但与无类函数相比,这会显着减慢该函数的速度(如结果所示(。因此,有没有另一种方法可以实现基于类的函数而不必调整向量的大小?为什么无类函数总是有正确的f
大小,而类函数没有?
在类求解函数中定义x_0
如下:
state_type x_0(dim, 1.);
并从类求解器的operator()
中删除resize
结果将如下所示:
GSL-error: 1.40024e-06, solved within 31397 ms
Boost-error: 1.74368e-08, solved within 152 ms
Boost_class-error: 1.74368e-06, solved within 1985 ms
仍然慢很多,但这是由于结构创建的开销?开销,如果重用类,则可以避免开销。
相关文章:
- 在C++中调整向量中的索引
- 将 boost::odeint 与向量类一起使用,而无需调整向量的大小
- 如何按另一个向量的方向调整一个向量?
- C++初始化列表中的向量集大小或调整大小
- std::向量迭代器和调整大小/保留的奇怪/有趣行为
- 调整向量大小并检索值,这是否正确或在任何情况下都可能导致段错误?
- 如何调整作为结构或类成员的向量的大小?
- 调整STD ::向量的大小是否可以降低其能力
- 在具有容量/调整大小的类中初始化向量
- STD ::向量如何调整其内部缓冲区大小
- 将添加的列初始化为可调整大小的 2D 向量
- C++中向量容器的大小和调整大小函数
- 如何判断 std::vector 是否调整了自身大小,以及如何解释指向向量内值的指针不再有效
- 调整OpenCV Mat向量向量的大小时出乎意料的结果
- C++有没有一种方法可以将向量调整为特定的字符
- 如何调整结构向量的大小
- 犰狳与 Boost Odeint 冲突:Odeint 在集成期间将状态向量调整为零
- 调整C 的2D向量大小,并填充0
- BAD_ALLOC错误实现向量调整大小函数时
- 如何在使用向量调整大小时调用非默认构造函数