2017-07-08 95 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本身。 –