2017-07-08 99 views
0

我打算解决几个矩阵微分方程,形式为d/dt (X) = F(X),其中X是一个大的复数矩阵,F表示它的一些函数。我试图用Boost的odeintstate_type作为Armadillo的cx_mat。但它会为受控步进器类型生成编译错误。我的示例代码如下犰狳的cx_mat和Boost的odeint编译错误

#include <armadillo> 
#include <iostream> 
#include <boost/numeric/odeint.hpp> 

using namespace std; 
using namespace arma; 
using namespace boost::numeric::odeint; 

using state_type = arma::cx_mat; 


void ode(const state_type& X, state_type& derr, double) { 
    derr = X; // sample derivative, can be anything else 
} 

// define resizable and norm_inf 
namespace boost { namespace numeric { namespace odeint { 

template <> 
struct is_resizeable<arma::cx_mat> { 
    typedef boost::true_type type; 
    const static bool value = type::value; 
}; 

template <> 
struct same_size_impl<arma::cx_mat, arma::cx_mat> { 
    static bool same_size(const arma::cx_mat& x, const arma::cx_mat& y) 
    { 
     return arma::size(x) == arma::size(y); 
    } 
}; 

template<> 
struct resize_impl<arma::cx_mat, arma::cx_mat> { 
     static void resize(arma::cx_mat& v1, const arma::cx_mat& v2) { 
     v1.resize(arma::size(v2)); 
    } 
}; 


template<> 
struct vector_space_norm_inf<state_type> { 
    typedef double result_type; 
    result_type operator()(const state_type& p) const 
    { 
     return arma::norm(p, "inf"); 
    } 
}; 



} } } // namespace boost::numeric::odeint 





using stepper = runge_kutta_dopri5<state_type, double, state_type, double, vector_space_algebra>; 

int main() { 

    cx_mat A = randu<cx_mat>(4, 4); 


    integrate_adaptive(make_controlled<stepper>(1E-10, 1E-10), ode, A, 0.0 , 10.0, 0.1); 
} 

此代码提供了以下编译错误:

/usr/include/armadillo_bits/Mat_meat.hpp:5153:3: error: static assertion failed: error: incorrect or unsupported type 
    arma_type_check((is_same_type< eT, typename T1::elem_type >::no)); 

什么我可以理解,犰狳不支持复制一个真正的矩阵(mat)成一个复杂的问题(cx_mat ),像

mat Z = something; 
cx_mat Y = Z; // ERROR 

其中在odeint某处发生。现在我通过整个矩阵复制到std::vector<std::complex<double> >克服这一放入功能ode,然后在函数内部重新复制整个std::vector<std::complex<double> >cx_mat,计算F(X),然后将其复制到std::vector<std::complex<double> >和回报。显然,这是非常缓慢和低效的。

任何简单的解决方案?如果可能的话,我可能想转向Eigen,如果有帮助的话。

回答

0

是的,复杂值状态类型不适合自适应步进器,因为odeint不区分state_type和error_type。即它也使用state_type来存储错误,但这对于复数值的state_types不起作用,而错误应该是双值矩阵。我不确定Eigen在这方面的表现如何,但它可能值得一试。随意提交https://github.com/headmyshoulder/odeint-v2的票,但这将是一个相当大的变化...

+0

Eigen没有帮助。我认为问题在于odeint本身。 –