本文へスキップ
モンテカルロシミュレーション
6.ランダムウォークのシミュレーション

 液体や気体中に浮遊する微粒子は、不規則な運動をくりかえしています。この現象は、1827年に植物学者ブラウンによって、水に浮く花粉の不規則な運動 として初めて発見されました。この現象は、酔っぱらいが、ちどり足で歩くのとにていることから、ランダムウォークあるいは酔歩の問題と呼ばれます。

1次元ランダムウォーク
 まずはじめに、簡単な1次元のランダムウォークを考えてみます。図6.1.1のように1本の道の上で、酔っぱらいが、ふらふらといったりきたりしているとします。

図6.1.1  1次元ランダムウォーク



酔っぱらいが、ある時刻で x の場所にいたとして、次の時刻に前へ行くか、後ろへ行くかは、まったくランダムで、予測が着きません。このような場合の、 x について、シミュレーションしてみましょう。


プログラムのポイント
 プログラムs0610.cは、1次元ランダムウォークのシミュレーション・プログラムで、各クロックごとに、行26において確率1/2で、前へ進むか、後ろに進むかをきめています。

 その結果を図6.1.2に示します。横軸時刻 t 、縦軸 x のグラフとして表現しています。このランダムウォークは、3.1節でのべたベルヌーイ試行を時刻 t ごとにおこなっているものです。
 またこれは、勝ったり負けたりのゲームと考えれば、持ち金の金額の変化をあらわしています。

図6.1.2  1次元ランダムウォークのシミュレーション結果

s0610.c   1次元ランダムウォーク
/* s0610.c 
 *        1次元ランダムウォーク 
 *         (C) H.Ishikawa 1994 2018
 */

#include "compati.h"
#include "window.h"

#define T_END   1000                 /* 終りの時刻 */
#define P        0.5                 /* 確率 */

int main(void)
{
  long    t;
  long    x = 0;
  
  /* グラフィクの準備 */
  opengraph();
  SetWindow(0,  0.0,-100.0,T_END,100.0,  SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
  Axis(0, "t", T_END/10, "x", 10);
  MoveTo(0, 0.0, 0.0);

  /* メイン */
  srand(1);
  for (t = 0; t < T_END; t ++) {
    if (P > RAND()) {
      x = x + 1;
    } else {
      x = x - 1;
    }
    LineTo(0, (double)(t), (double)(x));
  }
  getch();
  closegraph();
  return 0;
}



2次元ランダムウォーク
 平面上をブラウン運動する粒子を考えます。周囲の流体の分子から乱雑な衝撃を受け、ブラウン粒子はランダムな方向に運動をしますが、その方向は、 0 から 2π の間に一様に分布すると考えられます。衝撃を受けたブラウン粒子は、そのたびごとにさまざまな距離を進むと考えられますが、ここでは、図6.2.1のよう に一定の距離を進むものとします。このような、2次元のランダムウォークをシミュレーションしてみましょう。
図6.2.1 2次元ランダムウォーク



プログラムのポイント
 プログラムs0620.cでは、行32において0 から 2π の乱数を発生させ、行33,34のように x, y から、次の時点の x, y へ一歩進ませています。行35でブラウン粒子の軌跡をディスプレイします。 シミュレーション結果を図6.2.2に示します。


図6.2.2  2次元ランダムウォークのシミュレーション結果

s0620.c    2次元ランダムウォーク
/* s0620.c 
 *        2次元ランダムウォーク 
 *         (C) H.Ishikawa 1994 2018
 */

#include <math.h>                    /* sin() cos() */
#include "compati.h"
#include "window.h"

#define T_END    10000               /* 終りの時刻 */
#define PI       3.141593
#define X_MAX    80.0
#define Y_MAX    60.0

int main(void)
{
  double x,y;                       /* 粒子の座標 */
  double th;                        /* θ */
  int    t;
  
  /*グラフィクの準備*/
  opengraph();
  SetWindow(0,  -X_MAX, -Y_MAX, X_MAX, Y_MAX, 0,HIGHT, WIDTH,0);
  Axis(0, "", 10,"", 10);
  MoveTo(0, 0, 0);

  /*メイン*/
  srand(3);
  x = 0.0;
  y = 0.0;
  for (t = 1; t <= T_END; t ++) {
    th = 2.0 * PI * RAND();
    x = x + cos(th);
    y = y + sin(th);
    LineTo(0, x, y);
  }
  getch();
  closegraph();
  return 0;
}



拡散過程 (文献(10)参照)

 たくさんの粒子を一ケ所に集めて、ブラウン運動をさせると、粒子はどんどん広がって、ついには、一様にひろがってしまうでしょう。このような過程を拡散過程といいます。本節では、この拡散過程をシミュレーションします。
 ところで、単位時間に1度ずつ、衝撃を受けるブラウン運動で、 t=0 で原点にいた粒子が、 t=n にはどこにいる確率が高いでしょうか。図6.3.1のように、一回の衝撃で動く距離を1とすると、ブラウン粒子が、 n 回目の衝撃をうけると、原点からの距離 dn は、つぎのように求められます。
・・・・式6.3.1

となり、 n を大きくしていくと、平方根の中の2項目は n にくらべて無視できるようになります。したがって
・・・・・ 式6.3.2
となり、原点からの平均距離は衝撃の回収がふえると、その回数の平方根に比例して遠くなることがわかります。
 たくさんの粒子を一ケ所に集めて、ブラウン運動をさせると、粒子はどんどん広がって、ついには一様に分布してしまいます。このことを、実際にシミュレーションで確かめてみましょう。
図6.3.1 粒子の位置


プログラムのポイント
 プログラムs0630.cにおいては、2つのウインドウを用います。ウインドウ0では、粒子が運動する平面、ウインドウ1は、原点からの粒子の平均距離を表示することとします。
 時刻 0 で、2000個の粒子は原点にあるとし、クロックを1進めるたびに、前節と同じように、 0〜2π の乱数を発生させ、ランダムウォークをさせます。2000個の粒子に対しこのことを実行し、原点からの粒子の平均距離 d を計算したあと、時刻を進めていきます。


実行結果
 このプログラムを実行してみましょう。図6.3.2のように、はじめ原点にあった粒子は、あたかも生き物のように動きまわり、全体はしだいに拡散してい く様子が表示されます。また、たしかに原点からの平均距離は、ウインドウ1に示すように、時刻の平方根に比例していることも示されます。


図6.3.2 拡散のシミュレーション結果(途中の図)
時刻=10
時刻=100
時刻=500

  s0630.c   拡散のシミュレーション(長時間CPU占有に注意)
/* s0630.c 
 *        拡散のシミュレーション 
 *         (C) H.Ishikawa 1994 2018
 */

#include <math.h>                           /* sin(), cos(), sqrt() */
#include "compati.h"
#include "window.h"

#define PI       3.141593
#define PARTICLE 2000                       /* 粒子の数 */
#define MAX      50.0
#define T_END    500

double    x[PARTICLE] = {0.0};              /* すべての粒子は原点にある */
double    y[PARTICLE] = {0.0};
int       i;                                /* 粒子のカウンタ */
int       t = 0;
double    th;                               /* θ */
double    d1;                               /* 原点からの距離の和 */
double    d;                                /* 原点からの距離の平均値 */

int main(void)
{
  /*グラフィックの準備*/
  opengraph();
  SetWindow(0, -MAX,-MAX,MAX,MAX,  SPACE,HIGHT-SPACE,HIGHT-SPACE,SPACE);
  Axis(0, "", MAX, "", MAX);
  SetWindow(1,  0,0,T_END,MAX,  HIGHT,HIGHT/2,WIDTH-SPACE,SPACE);
  Axis(1, "", T_END / 5, "", MAX / 5);
  MoveTo(1, 0.0, 0.0);
  
  /*メイン*/
  while (t < T_END) {
    d1 = 0.0;
    for (i = 0; i < PARTICLE; i ++) {
      th = 2.0 * PI * RAND();
      PutPixel(0, x[i], y[i], 0);           /* 前の粒子の位置を消す */
      x[i] = x[i] + cos(th);
      y[i] = y[i] + sin(th);
      PutPixel(0, x[i], y[i], 7);           /* 次の粒子の位置を表示 */
      d1 = d1 + sqrt(x[i] * x[i] + y[i] * y[i]);
    }
    Axis(0,  "",MAX,"",MAX);                /* 消された軸を再度書く */
    d = d1 / PARTICLE;
    t = t + 1;
    LineTo(1,  t,d);
  }
  getch();
  closegraph();
  return 0;
}



周辺部では遅く動く
 前節のプログラムを長時間まわしていると、実は変なことが起こります。中心部にあった粒子が徐々に拡散していき、ついに MAX という枠をはみだしてしまいます。そこで、今回は枠まで粒子が到着したら、そこに壁があり、その壁ではねかえるようにします。
 また、前節では、粒子が1クロックで動く長さを一定としていましたが、中心部では速く、周辺部では遅く動くようにしてみましょう。どんなことがおこるでしょうか 。

プログラムのポイント
 プログラムs0640.cにおいては、1000個の粒子をブラウン運動させます。ただし、初期状態として一様に分布させることとし、また粒子の動き方は、1クロックで進む長さを l とすれば、
・・・・・ 式6.4.1
とします( MAX は壁までの距離、 xi, yi は粒子 i の位置)。すなわち、中心部で最も速く、
・・・・・ 式6.4.2
でうごき、中心部から最も遠い隅の
・・・・・ 式6.4.3
の時、停止状態となるものとします。
 また、次のクロックにおける位置を計算した結果、その値が MAX をこえた場合、図6.4.1のように、壁のところで、完全反射するものとします。すなわち
  if (x[i] >= MAX) { x[i] = MAX - (x[i] - MAX) ; }

のように MAX からはみだした分、ひき返すように計算します。 y 座標も同様です。

図6.4.1 粒子の動き方

実行結果
 このプログラムを実行してみましょう。部屋のほこりは、中心部では風が吹いたり、人の動きがあるため、速く動きまわり、次第に風の弱い部屋の隅に集まってきます。このような様子が図6.4.2のように、シミュレーションによって確かめることができます。


図6.4.2 ほこりはすみにたまるシミュレーション結果
開始時には均一に分布していたホコリは
時間がたつと下のように隅にたまっていく

開始時

時刻=5000

  s0640.c    ホコリのシミュレーション
/* s0640.c 
 *        ホコリのシミュレーション 
 *         (C) H.Ishikawa 1994 2018
 */

#include <math.h>                           /* sin(),  cos(), sqrt() */
#include "compati.h"
#include "window.h"

#define PI       3.141593
#define PARTICLE 1000                       /* 粒子の数 */
#define MAX      50.0
#define T_END    50000

double    x[PARTICLE];
double    y[PARTICLE];
int       i;                                /* 粒子のカウンタ */
long      t = 0;
double    th;
double    l;

int main(void)
{
  /*グラフィックの準備*/
  opengraph();
  SetWindow(0, -MAX,-MAX,MAX,MAX,  SPACE,HIGHT-SPACE,HIGHT-SPACE,SPACE);
  Axis(0, "", MAX, "", MAX);
  
  for (i = 0; i < PARTICLE; i ++) {
    x[i] = MAX * (2.0 * RAND() - 1.0);      /* 始めは一様に分布 */
    y[i] = MAX * (2.0 * RAND() - 1.0);
  }
  
  /*メイン*/
  while (t < T_END) {
    for (i = 0; i < PARTICLE; i ++) {
      PutPixel(0, x[i], y[i], 0);
      th = 2.0 * PI * RAND();
      l = (sqrt(2.0) * MAX - sqrt(x[i] * x[i] + y[i] * y[i])) / MAX;
      x[i] = x[i] + l * cos(th);
      y[i] = y[i] + l * sin(th);
      if (x[i] >=  MAX) { x[i] =  2*MAX - x[i];}   /* 反射 右の壁 */
      if (x[i] <= -MAX) { x[i] = -2*MAX - x[i];}   /*       左の壁 */
      if (y[i] >=  MAX) { y[i] =  2*MAX - y[i];}   /*       上の壁 */
      if (y[i] <= -MAX) { y[i] = -2*MAX - y[i];}   /*       下の壁 */
      PutPixel(0, x[i], y[i], 7);
    }
    Axis(0, "", MAX, "", MAX);
    t = t + 1;
  }
  getch();
  closegraph();
  return 0;
}



. プログラムs0620.cでは、歩幅は一定であったが、歩く度に歩幅がランダムに変わる場合をシミュレーションしなさい。(時間間隔一定、歩幅不定のランダムウォーク)

  (解答) 
  a0610.c  6章1.解答   2次元ランダムウォーク (歩幅がランダムに変わる場合)
/* a0610.c   6章1.解答
 *        2次元ランダムウォーク (歩幅がランダムに変わる場合)
 *         (C) H.Ishikawa 1994 2018
 */

#include <math.h>                    /* sin() cos() */
#include "compati.h"
#include "window.h"

#define T_END    10000               /* 終りの時刻 */
#define PI       3.141593
#define X_MAX    80.0
#define Y_MAX    60.0

int main(void)
{
  double x,y;                       /* 粒子の座標 */
  double th;                        /* θ */
  double l;                         /* 歩幅 */
  int    t;
  
  /*グラフィクの準備*/
  opengraph();
  SetWindow(0,  -X_MAX, -Y_MAX, X_MAX, Y_MAX, 0,HIGHT, WIDTH,0);
  Axis(0, "", 10,"", 10);
  MoveTo(0, 0, 0);

  /*メイン*/
  srand(31);
  x = 0.0;
  y = 0.0;
  for (t = 1; t <= T_END; t ++) {
    th = 2.0 * PI * RAND();
    l = -1.0 * log(1 - RAND());     /* 平均1.0の指数乱数 */
    x = x + l * cos(th);
    y = y + l * sin(th);
    LineTo(0, x, y);
  }
  getch();
  closegraph();
  return 0;
}


.  p が小さい場合のベルヌーイ試行をおこない、事象が起ったときに再度確率1/2で、右に行くか左に行くかきめるとします。このときの、軌跡を図示するプログラムを作りなさい。(歩幅一定、あるく時間が不定のランダムウォーク)

  (解答)  
  a0620.c 6章2.解答   1次元ランダムウォーク (歩幅一定,時間不定)
/* a0620.c   6章2.解答
 *        1次元ランダムウォーク (歩幅一定,時間不定)
 *         (C) H.Ishikawa 1994 2018
 */

#include "compati.h"
#include "window.h"

#define T_END   1000                /* 終りの時刻 */
#define P0      0.1                 /* 確率0 */
#define P1      0.5                 /* 確率1 */

int main(void)
{
  long    t;
  long    x = 0;
  
  /* グラフィクの準備 */
  opengraph();
  SetWindow(0,  0.0,-100.0,T_END,100.0,  SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
  Axis(0, "t", T_END/10, "x", 10);
  MoveTo(0, 0.0, 0.0);

  /* メイン */
  srand(3);
  for (t = 0; t < T_END; t ++) {
    if (P0 > RAND()) {              /* 小さい確率のベルヌーイ試行 */
      if (P1 > RAND()) {
        x = x + 1;
      } else {
        x = x - 1;
      }
    }
    LineTo(0, (double)(t), (double)(x));
  }
  getch();
  closegraph();
  return 0;
}