在 boost/odeint 中使用多个特征矩阵作为状态类型

Question

我正在研究一个物理问题，为此我必须根据 ODE 演化参数。他们必须不时地进行操作，以便我想要一种可以与对角化等例程一起使用的数据类型，...因此，我实现了一个 class 和 eigen::Matrix 作为成员和想要执行与 odeint 的集成。对于单个 eigen::matrix 这很好用。我做了一个最小的例子：

#include <iostream>
#include <fstream>
#include <cmath>
#include <string>

#include <Eigen/Core> 
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>

// define vector_space_algebra for Eigen::Matrix
namespace boost::numeric::odeint {
  template<typename B,int S1,int S2,int O, int M1, int M2>
  struct algebra_dispatcher< Eigen::Matrix<B,S1,S2,O,M1,M2> >{
    typedef vector_space_algebra algebra_type;
  };
}

// define abs() for Eigen::Matrix
namespace Eigen {
  template<typename D, int Rows, int Cols>
  Matrix<D, Rows, Cols> abs(Matrix<D, Rows, Cols> const& m) {
    return m.cwiseAbs();
  }
}

typedef Eigen::Matrix<double, 3,3> mat;
using namespace Eigen;
using namespace std;

class state {
 public:

  // state components
  Eigen::Matrix<double, 3,3> M1, M2;

  // constructors
  state() : M1(), M2() {};      // constructors
  state(mat M1in, mat M2in) : M1(M1in), M2(M2in) {};


  // in place addition and multiplication
  state operator+=(const state & X){
    M1 += X.M1; M2 += X.M2;
    return *this;
  }

  state operator*=(const double a){
    M1 *= a; M2 *= a;
    return *this;
  }

  // ODE
  void operator() ( const state & X , state & dX, const double ){
    dX.M1 = X.M1*X.M2.adjoint()*X.M2;
    dX.M2 = X.M2*X.M1.adjoint()*X.M1;
  }
};


// vector space operations

state operator+( const state &lhs , const state &rhs ){
  return state( lhs.M1+rhs.M1 ,lhs.M2+rhs.M2);
}

state operator*( const state &lhs , const double &rhs ){
  return state( lhs.M1*rhs ,lhs.M2*rhs);
}

state operator*( const double &lhs , const state &rhs ){
  return state( lhs*rhs.M1 ,lhs*rhs.M2);
}

state operator/( const state &lhs , const state &rhs ){
  return state( lhs.M1.cwiseQuotient(rhs.M1), lhs.M2.cwiseQuotient(rhs.M2) );
}
state abs( const state &X ){
  return state( abs(X.M1) , abs(X.M2) );
}


// lp infinity norm
namespace boost::numeric::odeint {
  template<>
  struct vector_space_norm_inf< state > {
    typedef double result_type;
    double operator()( const state &X ) const {
      return max( X.M1.lpNorm<Infinity>() , X.M2.lpNorm<Infinity>() );
    }
  };
}


//write to std output
void write( state &x , const double t ){
  cout << t << "\t" << x.M1 << "\t" << x.M2 << "\n";
}


//
// int main
//
int main(int argc, char* argv[]){

  // set values
  mat M1, M2;
  double t_end = 1;
  double t_start = 10;

  M1 << 0.1,0,0, 0,0.2,0.1 ,0.2,0,0.3;
  M2 << 0.5,0,0, 0,0.6,0, 0,0,0.7;

  state values(M1,M2);

  using namespace boost::numeric::odeint;

  // type definition for numerical integration
  typedef runge_kutta_dopri5< state , double, state , double,     vector_space_algebra > stepper;

  // integration
  int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10     ) , state() , values , t_start , t_end , 0.01);

  //output
  write(values,t_end);
  return(0);
}

基本上，这是取自 here

的示例

当我注释掉以“int steps”开头的行时，mac 上的 g++（我知道，那是不同的东西，错误更易读......）编译没有错误。否则，我得到

In file included from minimal.cpp:7:
In file included from /usr/local/include/boost/numeric/odeint.hpp:25:
In file included from /usr/local/include/boost/numeric/odeint/util/ublas_wrapper.hpp:30:
/usr/local/include/boost/numeric/odeint/algebra/default_operations.hpp:443:76: error: 
      invalid operands to binary expression ('double' and 'state')
  ...m_eps_abs + m_eps_rel * ( m_a_x * abs( get_unit_value( t1 ) ) + m_a_dxdt * abs( get_unit_value( t2 ) ) )...

因为它没有在那个文件中声明，所以我不明白函数 get_unit_value() 做了什么或想要我做什么。它似乎与误差估计有关，或者至少与执行集成中的某个步骤有关。我该如何解决？

Answer 1

问题是您没有在 double 和 state 之间定义 operator+。添加以下内容，代码至少应该可以编译。

state operator+( const state &lhs , double rhs ){
  return state( lhs.M1+rhs ,lhs.M2+rhs);
}

state operator+( double lhs , const state &rhs ){
  return state( lhs+rhs.M1 ,lhs+rhs.M2);
}

在 boost/odeint 中使用多个特征矩阵作为状态类型

using several eigen matrices as statetypes in boost/odeint

c++

eigen3

odeint