boost/odeintのstatetypeとしていくつかの固有行列を使用

Question

私はODEに従ってパラメータを展開しなければならない物理問題に取り組んでいます。対角化などのルーチンで使用できるデータ型を持つように、時にはそれらを操作する必要があります。したがって、メンバーとして固有クラス::行列を実装し、統合を実行したいodeintと一緒に。単一の固有値::行列に対して、これはうまくいきました。 、私が知っている（Macの場合はG ++、基本的に、これは私が「の手順をint型」で始まる行をコメントアウトするとhere

から取った例である

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

#include <eigen/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 { namespace numeric { namespace 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 { namespace numeric { namespace 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); 
} 

別の何かのthats：私は、最小限の例を作りましたどのようなエラーが読みやすく...）エラーなしでコンパイルされます。それ以外の場合は、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

operator+をdoubleとstateの間に定義していないという問題があります。次のコードを追加し、少なくともコンパイルする必要があります。

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); 
}