Eigen::VectorXd 和 Boost::Odeint，不起作用

Eigen::VectorXd and Boost::Odeint, not working

本文关键字：不起作用 Odeint VectorXd Eigen Boost      更新时间：2023-10-16

我正在使用Eigen：：VectorXd作为boost：：odeint的state_type进行测试。我使用此代码：

#include <Eigen/Eigen>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>
#include <iostream>
#include <vector>
template<class T>
struct push_back_state_and_time
{
std::vector<T>& m_states;
std::vector< double >& m_times;
push_back_state_and_time( std::vector<T> &states ,std::vector<double> &times )
: m_states(states) , m_times(times) { }
void operator()(const T &x ,double t )
{
m_states.push_back(x);
m_times.push_back(t);
}
};
template<class T>
struct write_state
{
void operator() ( const T &x, const double t ) {
std::cout << t << "t";
for(size_t i = 0; i < x.size(); ++i)
std::cout << x[i] << "t";
std::cout << std::endl;
};
};
template<class T>
class odeClass {
private:
double Ka, Kel, Vd;
public:
odeClass(double ka, double kel, double vd) : Ka(ka), Kel(kel), Vd(vd) {};
void operator() (const T &x, T &dxdt, const double t) {
dxdt[0] = - Ka * x[0];
dxdt[1] = Ka * x[0] - Kel * x[1];
};
};
void testODE_Eigen() {
double Ka = 0.195, Vd = 13.8, Kel = 0.79 / Vd;
Eigen::VectorXd x(2);
x << 40 / Vd, 0.0;
odeClass<Eigen::VectorXd> myClass(Ka, Kel, Vd);
boost::numeric::odeint::runge_kutta4<Eigen::VectorXd, double, Eigen::VectorXd, double, boost::numeric::odeint::vector_space_algebra> stepper;
size_t steps = integrate_adaptive( stepper, myClass, x ,0.0 ,100.0 ,0.5 ,write_state<Eigen::VectorXd>() );
}
void testODE_Vector() {
double Ka = 0.195, Vd = 13.8, Kel = 0.79 / Vd;
std::vector<double> x = { 40 / Vd, 0.0 };
odeClass<std::vector<double>> myClass(Ka, Kel, Vd);
boost::numeric::odeint::runge_kutta4<std::vector<double>> stepper;
size_t steps = integrate_adaptive( stepper, myClass, x ,0.0 ,100.0 ,0.5 ,write_state<std::vector<double>>() );
}
int main()
{
testODE_Eigen();
return 0;
}

运行时函数testODE_Vector();, it works perfectly, but when runningtestODE_Eigen((;' 我得到了第一个集成步骤，一个断言停止："断言失败：索引>= 0 &&索引<size((，文件>

有什么线索吗？

谢谢。

看起来您正在使用的特征库内的断言失败了。对于像Assertion failed: index >= 0 && index < size()这样的消息，库可能正在尝试在内部迭代向量，并在尝试访问它之前检查该向量是否有效。我会检查您传入的对象并确保它们有效。

看起来两个功能的主要区别之一是testODE_Vector()和testODE_Eigen()是您创建该x的方式。我不确定这行是什么x << 40 / Vd, 0.0;的目的是这样做，但我会从那里开始，并在x的值传递到integrate_adaptive之前

验证它是否正确

我的回答有点晚了，但万一其他人遇到这个问题，这是我发现的。

问题似乎是 OdeInt 无法正确处理特征向量和矩阵的动态大小。因此，在创建 dxdt 时，它会创建一个空的动态矩阵或向量。这会导致运算符重载中出现问题，即您尝试访问不包含任何 dxdt 的元素。

我发现的一个快速解决方法是使用resize()函数(或conservativeResize()(来确保 dxdt 具有适当的大小：

void operator() (const T &x, T &dxdt, const double t) {
dxdt.resize(x.size())
dxdt[0] = - Ka * x[0];
dxdt[1] = Ka * x[0] - Kel * x[1];
};

请注意，如果要使用矩阵而不是向量，则必须使用x.rows()和x.cols()而不是x.size()。

相关文章：