2.3.2 プログラムの解読

ball-bound-glsc.cpp

   1 /*
   2  * ball-bound-glsc.cpp --- Eigen を使ってはずむボールのシミュレーション
   3  *   http://nalab.mind.meiji.ac.jp/~mk/program/misc/ball-bound-glsc.cpp
   4  *   cglsc++ とか作れば良いのだけど
   5  *   c++ -I ~/include -I/opt/X11/include ball-bound-glsc.cpp -L ~/lib -lglscd -L/opt/X11/lib -lX11
   6  *   c++ -I/opt/X11/include ball-bound-glsc.cpp -lglscd -L/opt/X11/lib -lX11
   7  *
   8  *   ./a.out
   9  */
  10 
  11 #include <iostream>
  12 #include <math.h>
  13 #include <Eigen/Dense>
  14 using namespace Eigen;
  15 
  16 // GLSCのヘッダーを読み込む
  17 #ifndef G_DOUBLE
  18 #define G_DOUBLE
  19 #endif
  20 extern "C" {
  21 #include <glsc.h>
  22 };
  23 
  24 double m, g, Gamma, e;
  25 
  26 VectorXd f(double t, VectorXd x)
  27 {
  28   VectorXd y(4);
  29   y(0) = x(2);
  30   y(1) = x(3);
  31   y(2) = - Gamma / m * x(2);
  32   y(3) = - g - Gamma / m * x(3);
  33   return y;
  34 }
  35 
  36 int main(void)
  37 {
  38   int n, N;
  39   double tau, Tmax, t,pi;
  40   VectorXd x(4),k1(4),k2(4),k3(4),k4(4);
  41 
  42   pi = 4 * atan(1.0);
  43   m = 100;
  44   g = 9.8;
  45   Gamma = 1.0;
  46   e = 1.0;
  47 
  48   Tmax = 20;
  49   N = 1000;
  50   tau = Tmax / N;
  51   x << 0,0,50*cos(pi*50/180),50*sin(pi*50/180);
  52 
  53  /* グラフィックス・ライブラリィ GLSC の呼び出し */
  54   g_init((char *) "Meta", 250.0, 160.0);
  55   g_device(G_BOTH);
  56   g_def_scale(0, 0.0, 600.0, 0.0, 100.0, 10.0, 10.0, 230.0, 140.0);
  57   g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
  58   g_sel_scale(0);
  59   g_sel_line(0);
  60   g_cls();
  61   g_move(x(0),x(1));
  62 
  63   for (n = 0; n < N; n++) {
  64     t = n * tau;
  65     k1 = tau * f(t, x);
  66     k2 = tau * f(t+tau/2, x+k1/2);
  67     k3 = tau * f(t+tau/2, x+k2/2);
  68     k4 = tau * f(t+tau, x+k3);
  69     x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
  70     if (x(1)<0) {
  71       x(1) = - x(1);
  72       x(3) = - x(3);
  73     }
  74     std::cout << x(0) << " " << x(1) << std::endl;
  75     g_plot(x(0),x(1));
  76   }
  77   g_sleep(-1.0);
  78   g_term();
  79   return 0;
  80 }

• 11行目 これは C++ のプログラムには(実際上)いつでも書く。
• 12行目 C の #include <math.h> に相当する。 代わりに #include <math.h> としても OK (だと思う).
• 13行目 これで Eigen の密行列関係の機能が使えるようになる。
• 14行目 プログラムで以下 Eigen:: を省略するため。同様に using namespace std; とすることも考えられるが、 このプログラムでは、74行目に std::cout, std::endl としているだけなので、しなかった。
• 20〜22行 glsc.h は C 言語向けに書かれたヘッダー・ファイルであり、 C言語のプログラムでは #include <glsc.h> とすればインクルード出来るが、 C++ のプログラムでインクルードするには、このようなお呪いを前後につける必要がある。
• 26行目 VectorXd は、 成分が double で、サイズが動的に選べるベクトルを表す (Eigenの機能)。
• 26〜34行 関数定義。f() は、 常微分方程式 $\dfrac{\D\bm{x}}{\D t}=\bm{f}(t,\bm{x})$ の右辺の $\bm{f}(t,\bm{x})$ を計算するための関数である。
• 51行目。これは

    x(0) = 0;
    x(1) = 0;
    x(2) = 50 * cos(pi*50/180);
    x(3) = 50 * sin(pi*50/180);
  

  としても良い。
• 65〜69行はいわゆる Runge-Kutta法である。 時間刻みを $\tau$ としたとき、

  $$\bm{x}_{n+1}=\bm{x}_n+\frac{1}{6} \left(\bm{k}_1+2\bm{k}_2+2\bm{k}_3+\bm{k}_4\right),$$

  $$\bm{k}_1=\tau f(t,\bm{x})$$

  $$\bm{k}_2=\tau f(t+\frac{\tau}{2},\bm{x}+\bm{k}_1/2)$$

  $$\bm{k}_3=\tau f(t+\frac{\tau}{2},\bm{x}+\bm{k}_2/2)$$

  $$\bm{k}_4=\tau f(t+\tau,\bm{x}+\bm{k}_3)$$

  
  
  というのが Runge-Kutta法の公式であるが、 それをほぼそのまま書けるところが、ベクトルを使える利点である。