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

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

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

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



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


プログラムのポイント
 プログラムS0610.javaは、1次元ランダムウォークのシミュレーション・プログラムで、各クロックごとに、(1)において確率1/2で、前へ進むか、後ろに進むかをきめています。その結果を図6.1.2に示します。横軸時刻 t 、縦軸 x のグラフとして表現しています。このランダムウォークは、3.1節でのべたベルヌーイ試行を時刻 t ごとにおこなっているものです。
 またこれは、勝ったり負けたりのゲームで、持ち金の金額の変化をあらわしています。

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

S0610.java 1次元ランダムウォーク  | Download | 実行 |
/* S0610.java
 *        1次元ランダムウォーク  
 *         (C) H.Ishikawa 2008 
 */

package simulation;

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;


public class S0610 extends Applet implements ActionListener {
        Button button0;
        public void init() {
                button0 = new Button(" 再実行 ");
                add(button0);
                button0.addActionListener(this);
        }
        public void actionPerformed(ActionEvent e) {
                String label = e.getActionCommand();
                repaint();
        }

        public void paint(Graphics g) {
                Window w ;
                w = new Window();
                int SPACE = 30;
                int HIGHT = 400;
                int WIDTH = 640;
                
                long      t = 0;
                long      x = 0;
                long      T_END =   1000;    /* 終りの時刻 */
                double    P =     0.5;       /* 確率1/2 */

                
                /*グラフィックの準備*/
                w.setWindow(0,  0.0,-100.0,T_END,100.0,
                     SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
                w.axis(0, "t", T_END/10, "x", 10, g);
                
                w.moveTo(0, 0.0, 0.0, g);
  
                /*メイン*/
                for (t = 0; t < T_END; t ++) {
                        if (P > Math.random()) {   
                        //(1) 確率1/2で前へ進むか、後ろに進むか
                                x = x + 1;
                        } else {
                                x = x - 1;
                        }
                        g.setColor(Color.green);
                        w.lineTo(0, (double)(t), (double)(x), g);
                }
                stop();
        }
}


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



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


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



S0620.java 2次元ランダムウォーク  | Download | 実行 |
/* S0620.java
 *        2次元ランダムウォーク  
 *         (C) H.Ishikawa 2008 
 */

package simulation;

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;


public class S0620 extends Applet implements ActionListener {
	Button button0;
	public void init() {
		button0 = new Button(" 再実行 ");
		add(button0);
		button0.addActionListener(this);
	}
	public void actionPerformed(ActionEvent e) {
		String label = e.getActionCommand();
		repaint();
	}

	public void paint(Graphics g) {
		Window w ;
		w = new Window();
		int        SPACE = 30;
		int        HIGHT = 400;
		int        WIDTH = 640;
		
		double     X_MAX =   80.0;
		double     Y_MAX =   50.0;
		
		double x,y;                       /* 粒子の座標 */
		double th;                        /* θ */
		int    t;			  /*時刻*/
		long T_END =   10000;             /* 終りの時刻 */
		
		/*グラフィックの準備*/
		w.setWindow(0,  -X_MAX, -Y_MAX, X_MAX, Y_MAX,
                      SPACE,HIGHT-SPACE,WIDTH-SPACE,SPACE);
		w.axis(0, "", 10,"", 10, g);		
		w.moveTo(0, 0, 0, g);
  
		/*メイン*/
		x = 0.0;
		y = 0.0;
		for (t = 1; t <= T_END; t ++) {
			th = 2.0 * Math.PI * Math.random(); //(1)
		    x = x + Math.cos(th);               //(2)
		    y = y + Math.sin(th);               //(2)
		    g.setColor(Color.green);
		    w.lineTo(0, x, y, g);               //(3)
		}	
		stop();
	}
}


拡散過程 (文献(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.javaにおいては、2つのウインドウを用います。ウインドウ0では、粒子が運動する平面、ウインドウ1は、原点からの粒子の平均距離を表示することとします。
 時刻 0 で、1000個の粒子は原点にあるとし、クロックを1進めるたびに、前節と同じように、 0〜2π の乱数を発生させ、ランダムウォークをさせます。1000個の粒子に対しこのことを実行し、原点からの粒子の平均距離 d を計算したあと、時刻を進めていきます。


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


図6.3.2 拡散のシミュレーション結果(途中の図 長時間CPU占有注意)

S0630.java 拡散のシミュレーション  | Download | 実行 (長時間CPU占有注意)|
/* S0630.java
 *        拡散のシミュレーション 
 *         (C) H.Ishikawa 2008 
 */

package simulation;

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;


public class S0630 extends Applet implements ActionListener {
        Button button0;
        Button button1;
        int sw =0;
        public void init() {
                button0 = new Button(" 実行 ");
                button1 = new Button("クリア");
                add(button0);
                add(button1);
                button0.addActionListener(this);
                button1.addActionListener(this);
        }
        public void actionPerformed(ActionEvent e) {
                String label = e.getActionCommand();
                if (label.equals(" 実行 ")) sw = 1; else sw = 0;
                repaint();
        }

        public void paint(Graphics g) {
                Window w ;
                w = new Window();
                int SPACE = 30;
                int HIGHT = 400;
                int WIDTH = 640;
                int PARTICLE =  10000;                 /* 粒子の数 */
                double MAX  =   50.0;
                long T_END =   1000;

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

                /*グラフィックの準備*/
                w.setWindow(0, -MAX,-MAX,MAX,MAX,
            SPACE,HIGHT-SPACE,HIGHT-SPACE,SPACE);
                w.axis(0, "", MAX, "", MAX, g);
                w.setWindow(1,  0,0,T_END,MAX, 
            HIGHT,HIGHT/2,WIDTH-SPACE,SPACE);
                w.axis(1, "時間", T_END / 5, "原点からの平均距離", MAX / 5, g);
                w.moveTo(1, 0.0, 0.0, g);
  
                /*メイン*/
                if ( sw == 1) {
                        while (t < T_END) {
                                d1 = 0.0;
                                for (i = 0; i < PARTICLE; i ++) {
                                  th = 2.0 * Math.PI * Math.random();
                                  if (x[i] <=  MAX && x[i] >= -MAX &&
                                      y[i] <=  MAX && y[i] >= -MAX) {
                                          g.setColor(Color.white); 
                                          /* 前の粒子の位置を消す */
                                           w.putPixel(0, x[i], y[i], g);}          
                                   x[i] = x[i] + Math.cos(th);
                                   y[i] = y[i] + Math.sin(th);
                                   if (x[i] <=  MAX && x[i] >= -MAX && 
                                       y[i] <=  MAX && y[i] >= -MAX) {
                                          g.setColor(Color.blue);
                                          /* 次の粒子の位置を表示 */
                                          w.putPixel(0, x[i], y[i], g); }          
                                   d1 = d1 + Math.sqrt(x[i] * x[i] + y[i] * y[i]);
                                }
                                /* 消された軸を再度書く */
                                w.axis(0,  "",MAX,"",MAX, g);                
                                d = (double)(d1 / PARTICLE);
                                t = t + 1;
                                w.lineTo(1,  t, d, g);
                        }
                        stop();
                }
        }
}


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

プログラムのポイント
 プログラムS0640.javaにおいては、S0630.javaと同様に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 ほこりはすみにたまるシミュレーション結果(長時間CPU占有注意)

S0640.java ホコリのシミュレーション  | Download | 実行(長時間CPU占有注意) |
/* S0640.java
 *        ホコリのシミュレーション 
 *         (C) H.Ishikawa 2008 
 */

package simulation;

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;


public class S0640 extends Applet implements ActionListener {
        Button button0;
        Button button1;
        int sw =0;
        public void init() {
                button0 = new Button(" 実行 ");
                button1 = new Button("クリア");
                add(button0);
                add(button1);
                button0.addActionListener(this);
                button1.addActionListener(this);
        }
        public void actionPerformed(ActionEvent e) {
                String label = e.getActionCommand();
                if (label.equals(" 実行 ")) sw = 1; else sw = 0;
                repaint();
        }

        public void paint(Graphics g) {
                Window w ;
                w = new Window();
                int SPACE = 30;
                int HIGHT = 400;
                int WIDTH = 640;
                int PARTICLE =  1000;                       /* 粒子の数 */
                double MAX  =     50.0;
                long T_END =   10000;

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

                /*グラフィックの準備*/
                w.setWindow(0, -MAX,-MAX,MAX,MAX,
                     SPACE,HIGHT-SPACE,HIGHT-SPACE,SPACE);
                w.axis(0, "", MAX, "", MAX, g);
                
                /*開始*/
                if ( sw == 1) {
                        for (i = 0; i < PARTICLE; i ++) {
                                /* 始めは一様に分布 */
                                x[i] = MAX * (2.0 * Math.random() - 1.0);     
                                y[i] = MAX * (2.0 * Math.random() - 1.0);
                                g.setColor(Color.green);
                                w.putPixel(0, x[i], y[i], g);
                        }
  
                        /*メイン*/
                        while (t < T_END) {
                                for (i = 0; i < PARTICLE; i ++) {
                                        g.setColor(Color.white);
                                        w.putPixel(0, x[i], y[i], g);
                                        th = 2.0 * Math.PI * Math.random();
                                        l = (Math.sqrt(2.0) * MAX -
                                        Math.sqrt(x[i] * x[i] + y[i] * y[i])) / MAX;
                                        x[i] = x[i] + l * Math.cos(th);
                                        y[i] = y[i] + l * Math.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];}   
                                        g.setColor(Color.blue);
                                        w.putPixel(0, x[i], y[i], g);
                                }
                                w.axis(0, "", MAX, "", MAX ,g);
                                t = t + 1;
                        }
                        stop();
                }
        }
}


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

  (解答)  A0610.java  | 実行 |

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

  (解答)  A0620.java  | 実行 |