2.1 計算だけするプログラム

計算だけして、軌跡は gnuplot で、というプログラム

ball-bound.cpp

   1 /*
   2  * ball-bound.cpp --- Eigen を使ってはずむボールのシミュレーション
   3  *
   4  *   c++ ball-bound.cpp
   5  *     もしかしたら -I/usr/local/include のような指定が必要かも。
   6  *   ./a.out > ball-bound.dat
   7  *   gnuplot
   8  *   gnuplot> plot "ball-bound.dat"
   9  */
  10 
  11 #include <iostream>
  12 #include <cmath>
  13 #include <Eigen/Dense>
  14 using namespace Eigen;
  15 
  16 double m, g, Gamma, e;
  17 
  18 VectorXd f(double t, VectorXd x)
  19 {
  20   VectorXd y(4);
  21   y(0) = x(2);
  22   y(1) = x(3);
  23   y(2) = - Gamma / m * x(2);
  24   y(3) = - g - Gamma / m * x(3);
  25   return y;
  26 }
  27 
  28 int main(void)
  29 {
  30   int n, N;
  31   double tau, Tmax, t,pi;
  32   VectorXd x(4),k1(4),k2(4),k3(4),k4(4);
  33 
  34   pi = 4 * atan(1.0);
  35   m = 100;
  36   g = 9.8;
  37   Gamma = 1.0;
  38   e = 1.0;
  39 
  40   Tmax = 20;
  41   N = 1000;
  42   tau = Tmax / N;
  43   x << 0,0,50*cos(pi*50/180),50*sin(pi*50/180);
  44   for (n = 0; n < N; n++) {
  45     t = n * tau;
  46     k1 = tau * f(t, x);
  47     k2 = tau * f(t+tau/2, x+k1/2);
  48     k3 = tau * f(t+tau/2, x+k2/2);
  49     k4 = tau * f(t+tau, x+k3);
  50     x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
  51     if (x(1)<0) {
  52       x(1) = - x(1);
  53       x(3) = - x(3);
  54     }
  55     std::cout << x(0) << " " << x(1) << std::endl;
  56   }
  57   return 0;
  58 }

curl -O http://nalab.mind.meiji.ac.jp/~mk/ball-bound.cpp
head ball-bound.cpp
c++ ball-bound.cpp
./a.out > ball-bound.dat
gnuplot

(curl でファイルを入手。head で先頭 10 行を表示して中身を確認。 c++ でコンパイル。実行形式 a.out の出力を ball-bound.dat に保存。gnuplot を起動。)

gnuplot のプロンプトに gnuplot> に対して、 plot "ball-bound.dat" とすると軌跡が描かれる(はず)。 gnuplot を終了するには、quit とすれば良い。

gnuplot> plot "ball-bound.dat"
gnuplot> quit
$