PaleAle011235の備忘録

経緯

制御系のシミュレーションをしたいが, 自分でソルバを書くのは面倒.
でも環境の問題で新しいライブラリを入れるのも大変.
boostに入ってるじゃないか！

という経緯でで, Boost.Numeric.Odeintを試してみました.

対象とした問題

単にライブラリの使い方の練習をしたいだけなので, 問題はなんでも構わなかったのですが,

1. 一般解が簡単なので, 検証がしやすい.
2. 方程式が簡単
3. 単に馴染みがある

という理由で調和振動子の運動方程式を選択しました.

調和振動子はバネと質点からなる振動子です.
運動方程式は

$m\frac{d^{2}x}{dt^{2}} = -kx$

mは質点の質量. kはばね定数です.

数値解析ではアルゴリズムは $\frac{du}{dt} = f(u,t)$ という形をしている微分方程式を取り扱うため, その形式に変形.

$\frac{dx}{dt} = v \\\ \frac{dv}{dt} = -\frac{k}{m}x$

念のために, $u = (x, v)^{\rm T}$ という2次元の変数を定義して変形させておくと,

$\frac{du}{dt} = \begin{pmatrix} 0 & 1\\\ -\frac{k}{m} & 0 \end{pmatrix} u$

今回は, 計算を簡単にするため $m =1$ , $k = 1$ としています.

Boost.Numeric.Odeintの使い方

Odeintは右辺の関数を登録することにより微分方程式を定義し, 数値解を計算してくれる.
まずはOdeintを利用するために必要なヘッダファイル.

#include <boost/numeric/odeint.hpp>

次に, Odeintがいる名前空間.

using namespace boost::numeric::odeint;

いよいよ, 微分方程式の右辺を定義する.
ここで, 関数の戻り値と引数はかならず一致している必要がある.

//相空間は2次元
using state_type = array<double,2>; 
 
//調和振動子の状態方程式
#define K 1.0 //バネ定数
#define M 1.0 //質量
void ocirator_system(const state_type& x, state_type& dot_x, const double t)
{
    //d^2/dt^2 = dv/dt = -k/m*x
    //           dx/dt = v
    dot_x[0] = -K/M*x[1];
    dot_x[1] = x[0];
}

定義の関数をシステムに登録し, 計算を実行させる.
integrate_constでは1ステップごとに固定時間だけシミュレーションが進行する.
このシミュレーションでは $t=0$ から出発し, 0.001づつ時間が経過していく.
Odeintには陰的ルンゲクッタなどの手法も実装されており, 1ステップごとの時間経過を動的に切り替えることでステップ数を減らし, 実行時間を低減させることも可能となっている.

    //初期状態
    state_type state0= {1,0};// x0 = 1, v0 = 0
 
    //解析手法を選択:4次の陽的ルンゲクッタ
    runge_kutta4<state_type> rk;
 
    //数値解析を実行
    //ステップごとにステップ幅を変えない
    integrate_const(rk,              //ステップごとの手法
                        ocirator_system, //状態方程式
                        state0,          //初期値
                        0.0,             //初期時刻
                        M_PI*2*3,        //終了時刻
                        0.001,           //ステップ幅
                        [&](const state_type &x, const double t) //ステップ毎に実行される関数.
                        {
                                        //標準出力へ
                            cout << t << "," << x[0] << "," << x[1] << endl;
                        });

すべて合わせると, こんなかんじ.

#include "iostream"
#include "array"
#include <boost/numeric/odeint.hpp>
 
#define K 1.0 //バネ定数
#define M 1.0 //質量
 
using namespace std;
using namespace boost::numeric::odeint;
 
//相空間は2次元
using state_type = array<double,2>; 
 
//調和振動子の状態方程式
void ocirator_system(const state_type& x, state_type& dot_x, const double t)
{
    //d^2/dt^2 = dv/dt = -k/m*x
    //           dx/dt = v
    dot_x[0] = -K/M*x[1];
    dot_x[1] = x[0];
}
 
int main()
{
    //初期状態
    state_type state0= {1,0};// x0 = 1, v0 = 0
 
    //解析手法を選択:4次の陽的ルンゲクッタ
    runge_kutta4<state_type> rk;
 
    //数値解析を実行
    //ステップごとにステップ幅を変えない
    integrate_const(rk,              //ステップごとの手法
                        ocirator_system, //状態方程式
                        state0,          //初期値
                        0.0,             //初期時刻
                        M_PI*2*3,        //終了時刻
                        0.001,           //ステップ幅
                        [&](const state_type &x, const double t) //ステップ毎に実行される関数.
                        {
                                        //標準出力へ
                            cout << t << "," << x[0] << "," << x[1] << endl;
                        });
    return 0;
}