先日
えすだぶさんの高射砲の本が面白くて、どれくらい当てるの難しいんだろうかと思ったのでざっくり書いて見たけど、全然あたんないですね。高射砲の初速度1000 m/sで、グラフの左から右にフラフラしながら飛んでる飛行機を撃ち落としたい pic.twitter.com/lvEJQDHBxW— 鹿部 等 (@c_curve1870) 2019年3月16日
というような計算をしたときに、軌道を計算するためにscipyのodeintを使って計算したところ、思ったより簡単に、ちゃんとした公式の元で積分が行えるのでちょっと感動しました。
そういうわけで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に変換できないって書いてましたね。
そんなメモ書きでした。
