3.4 C++言語＆Eigenによるプログラム例

C++ は C言語と似ているので、C言語を知っている人には、 敷居が低いであろう (Mac で、C言語のコンパイラーを使えるようにしてある人は、 同時に C++ のコンパイラーも使えるようになっているはず)。 便利なライブラリィが豊富にある。 ベクトル・行列計算のためには、Eigen というライブラリィがおすすめ。

Eigenのインストールは簡単 (A を見よ。 5分程度で済むので、うまく行かなければ気軽に質問して下さい。)。

C++は、入出力が C 言語とは違うが、それについては、 「なんとか C++ を使う」 を見ると良い。

Eigen には、VectorXd というベクトルを扱うクラスがある。

#include <Eigen/Dense>
using namespace	Eigen;
...

  VectorXd(2) x;

とすると、2次元ベクトル x が定義できる。 ベクトル x の成分 x(0), x(1) が使えるようになる。 これだけだと配列と変わらないが、 ベクトル演算ができるのが大きな利点である。 ここでは2次元であるが、何次元でも

    k1 = tau * f(t, x);
    k2 = tau * f(t+tau/2, x+k1/2);
    k3 = tau * f(t+tau/2, x+k2/2);
    k4 = tau * f(t+tau, x+k3);
    x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;

により Runge-Kutta 法の計算ができる。

rk2ex1++.cpp

/*
 * rk2ex1++.cpp (Runge-Kutta method for van der Pol equation)
 *   c++ -I /usr/local/include rk2ex1++.cpp
 *   ./a.out > rk2ex1.data  (分割数Nを入力. 例えば1000)
 *   gnuplot で
 *    splot "rk2ex1.data" with l
 *    plot "rk2ex1.data" with l, "rk2ex1.data" using 1:3 with l
 *    plot "rk2ex1.data" using 2:3 with l
 */

#include <iostream>
#include <math.h>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

double mu = 1.0;

int main(void)
{
  int i, N;
  double t, dt, pi;
  VectorXd x(2),k1(2),k2(2),k3(2),k4(2),f(double,VectorXd),x0(2);
  double Tmax;
  // 最終時刻
  Tmax = 50.0;
  // 時間刻み
  cout << "# N: "; cin >> N;
  dt = Tmax / N;
  // 初期値
  x0 << 0.1,0.1;
  t = 0.0;
  x = x0;
  cout << "# t x" << endl;
  cout << fixed; // この後はC言語の %f の形式で出力
  cout << t << x(0) << x(1) << endl;
  // Runge-Kutta法
  for (i = 0; i < N; i++) {
    k1 = dt * f(t, x);
    k2 = dt * f(t+dt/2, x+k1/2);
    k3 = dt * f(t+dt/2, x+k2/2);
    k4 = dt * f(t+dt, x+k3);
    x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    t = (i + 1) * dt;
    cout << t << " " << x(0) << " " << x(1) << endl;
  }
}

VectorXd f(double t, VectorXd x)
{
  VectorXd y(2);
  y(0) = x(1);
  y(1) = - x(0) + mu * (1.0 - x(0) * x(0)) * x(1);
  return y;
}

コンパイル&実行

% c++ -O3 -I /usr/local/include -o rk2ex1 rk2ex1++.cpp
% ./rk2ex1 > rk2ex1++.data
1000
%

この後は、これまでと同様にして、gnuplot を用いて可視化ができる。

ところで、GLSC (付録 C を見よ) を用いると、 C++のプログラムからグラフィックスが使える。 以下のプログラム中で、 名前がg_で始まっている関数が、GLSCの関数である。 初期化や座標系の設定などはごちゃごちゃしているが、 グラフを描く部分は g_move(x,y), g_plot(x,y) の2命令しか使っていないので、 難しくはないだろう。

rk2ex1++glsc.cpp

/*
 * rk2ex1++glsc.cpp (Runge-Kutta method for van der Pol equation)
 *   c++ -I /usr/local/include -I/opt/X11/include -o rk2ex1++glsc \
 *       rk2ex1++glsc.cpp -L/usr/local/lib -lglscd -L/opt/X11/lib -lX11 -lm
 *   ./rk2ex1++glsc  (分割数Nを入力. 例えば1000)
 */

#include <iostream>
#include <math.h>

#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

extern "C" {
#include <glsc.h>
};

void g_dump(const char *fname, Display *display, Window wid)
{
  char command[256];
  sprintf(command, "import -window %lu %s", wid, fname);
  system(command);
}

double mu = 1.0;

int main(void)
{
  int i, N;
  double t, dt, pi;
  VectorXd x(2),k1(2),k2(2),k3(2),k4(2),f(double,VectorXd),x0(2);
  double Tmax;
  double a = -3.0, b = 3.0, c = -3.0, d = 3.0;
  double win_width, win_height, w_margin, h_margin;
  // 最終時刻
  Tmax = 50.0;
  // 時間刻み
  cout << "# N: "; cin >> N;
  dt = Tmax / N;
  // 初期値
  x0 << 0.1,0.1;
  t = 0.0;
  x = x0;
  cout << "# t x" << endl;
  cout << fixed; // この後はC言語の %f の形式で出力
  cout << t << " " << x(0) << " " << x(1) << endl;
  // GLSCの初期化、座標系の設定
  win_width = 150.0; win_height = 150.0; w_margin = 10.0; h_margin = 10.0;
  g_init((char *)"LOTKA_VOLTERRA",
         win_width  + 2 * w_margin, win_height + 2 * h_margin);
  g_device(G_BOTH);
  g_def_scale(0, a, b, c, d,
	      w_margin, h_margin, win_width, win_height);
  g_sel_scale(0);
  // 線種の定義
  g_def_line(0, G_BLACK, 2, G_LINE_SOLID);
  g_def_line(1, G_BLUE,  1, G_LINE_SOLID);
  /* 座標軸を描く */
  g_sel_line(0);
  g_move(a, 0.0);  g_plot(b, 0.0);
  g_move(0.0, c);  g_plot(0.0, d);
  // ここからグラフ描き
  g_sel_line(1);
  g_move(x(0), x(1));
  // Runge-Kutta法
  for (i = 0; i < N; i++) {
    k1 = dt * f(t, x);
    k2 = dt * f(t+dt/2, x+k1/2);
    k3 = dt * f(t+dt/2, x+k2/2);
    k4 = dt * f(t+dt, x+k3);
    x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    t = (i + 1) * dt;
    g_plot(x(0), x(1));
    cout << t << " " << x(0) << " " << x(1) << endl;
  }
  // クリックされるまで待つ
  g_sleep(-1.0);
  // 画像の保存
  g_dump("rk2ex1++glsc.png", g_get_display(), g_get_window());
  // GLSCの終了
  g_term();
  return 0;
}

VectorXd f(double t, VectorXd x)
{
  VectorXd y(2);
  y(0) = x(1);
  y(1) = - x(0) + mu * (1.0 - x(0) * x(0)) * x(1);
  return y;
}

このプログラムを実行すると、 画面にウィンドウが現れて、そこに解軌道が描かれる。 計算終了後にそのウィンドウをクリックすると、 画像をファイル (rk2ex1++glsc.png) に保存してからプログラムが終了する。

図 6: van der Pol 方程式の解軌道