2019年4月7日日曜日

自作コンテナクラスでboost odeintのcontrolled stepperを使うメモ書き

先日

というような計算をしたときに、軌道を計算するためにscipyのodeintを使って計算したところ、思ったより簡単に、ちゃんとした公式の元で積分が行えるのでちょっと感動しました。
で、C95で公開したプログラムもodeintで実装できると助かるなぁと思ったんですが、scipy.integrate.odeintは引数がリストのみなので、Anderson-Walkerモデルみたいに13個もある式を真面目にやろうとすれば絶対どこかでindexを間違える未来が見えます。

そういうわけでscipyのodeintを調べているとboost::odeintというものがあり、こちらは自作のコンテナクラスも使えるということを知りました。自作コンテナクラスが使えるのであれば、メンバーにstd::vectorを持ち、それらのindexに対してsetter, getterを作れば、計算部分を書く時は構造体っぽく使えるので良さそうだな、と思いました(自作のPoint classも使えるんですが、operatorを何個も手打ちしたくなかったので今回はstd::vectorをメンバに持つクラスを作りました)。
ここまでの段階ではC++で常微分方程式:boost::odeint入門およびboost odeintのドキュメントを多く参考にしました。

さて、せっかくodeintのようなパッケージが使えるんですから、Controlled stepperを使って可変時間刻みの元で計算をしたいわけですが、

void lorenz(const state_type &x, state_type &dxdt, const double t)
{
  const double sigma(10.0);
  const double R(28.0);
  const double b(8.0 / 3.0);

  dxdt[0] = sigma * (x[1] - x[0]);
  dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
  dxdt[2] = -b * x[2] + x[0] * x[1];
};
using namespace boost::numeric::odeint;

int main()
{
  state_type x;//git hubにあるサンプルコードに従って作った自作のコンテナクラス
  x[0] = 5.0;
  x[1] = 10.0;
  x[2] = 10.0;
  state_type xold = x;

  BOOST_STATIC_ASSERT(is_resizeable::value == true);

  using base_stepper_type = runge_kutta_dopri5;
  auto Stepper = make_controlled(1e-6, 1e-6, base_stepper_type());

  double t = 0.;
  double dt_init = 1e-2;
  double dt = dt_init;
  double dt_max = 1.;
  double tmax = 10.;
  boost::numeric::odeint::controlled_step_result res;
  while (t < tmax)
  {

    xold = x;
    res = Stepper.try_step(lorenz, x, t, dt);
    if (res)
    {
      res = Stepper.try_step(lorenz, x, t, dt);
    }
    else
    {
    }
    if (dt > dt_max)
    {
      dt = dt_max;
    }
    std::cout << x.x() << x.y() << x.z() << std::endl;
  };
}

みたいなことをやると、
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:89:40: error: cannot convert 'boost::numeric::odeint::norm_result_type<AwState, void>::type {aka AwState}' to 'boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}' in return
         return algebra.norm_inf( x_err );
以下、大量のエラーを吐いて爆死します。で、これを回避したい、というのが今回の話題。
C++全然わからないマンなりにいろいろ試した結果、git hubにあるサンプルコード

template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;

public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;
 
public:
    my_vector( const size_t N )
...


template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;
public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;
    typedef vector::value_type value_type; //ここを追加

public:
    my_vector( const size_t N )
...

とすればなんの問題もなく通るようになりました。
cannot convert 'boost::numeric::odeint::norm_result_type<AwState, void>::type {aka AwState}' to 'boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}' in return
         return algebra.norm_inf( x_err );
よく見たら、typeをvalue_typeに変換できないって書いてましたね。
そんなメモ書きでした。

0 件のコメント:

コメントを投稿