 |
液体や気体中に浮遊する微粒子は、不規則な運動をくりかえしています。この現象は、１８２７年に植物学者ブラウンによって、水に浮く花粉の不規則な運動
として初めて発見されました。この現象は、酔っぱらいが、ちどり足で歩くのとにていることから、ランダムウォークあるいは酔歩の問題と呼ばれます。
１次元ランダムウォーク
まずはじめに、簡単な１次元のランダムウォークを考えてみます。図６．１．１のように１本の道の上で、酔っぱらいが、ふらふらといったりきたりしているとします。
| 図６．１．１ １次元ランダムウォーク |
 |
酔っぱらいが、ある時刻で x の場所にいたとして、次の時刻に前へ行くか、後ろへ行くかは、まったくランダムで、予測が着きません。このような場合の、 x について、シミュレーションしてみましょう。
プログラムのポイント
プログラムs0610.cは、１次元ランダムウォークのシミュレーション・プログラムで、各クロックごとに、行26において確率１／２で、前へ進むか、後ろに進むかをきめています。
その結果を図６．１．２に示します。横軸時刻 t 、縦軸 x のグラフとして表現しています。このランダムウォークは、３．１節でのべたベルヌーイ試行を時刻 t ごとにおこなっているものです。
またこれは、勝ったり負けたりのゲームと考えれば、持ち金の金額の変化をあらわしています。
| 図６．１．２ １次元ランダムウォークのシミュレーション結果 |
 |
| s0610.c １次元ランダムウォーク |
/* s0610.c
* １次元ランダムウォーク
* (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;
}
|
２次元ランダムウォーク
平面上をブラウン運動する粒子を考えます。周囲の流体の分子から乱雑な衝撃を受け、ブラウン粒子はランダムな方向に運動をしますが、その方向は、 0
から 2π
の間に一様に分布すると考えられます。衝撃を受けたブラウン粒子は、そのたびごとにさまざまな距離を進むと考えられますが、ここでは、図６．２．１のよう
に一定の距離を進むものとします。このような、２次元のランダムウォークをシミュレーションしてみましょう。
| 図６．２．１ ２次元ランダムウォーク |
 |
プログラムのポイント
プログラムs0620.cでは、行32において0 から 2π の乱数を発生させ、行33，34のように x, y から、次の時点の x, y へ一歩進ませています。行35でブラウン粒子の軌跡をディスプレイします。 シミュレーション結果を図６．２．２に示します。
| 図６．２．２ ２次元ランダムウォークのシミュレーション結果 |
 |
| s0620.c ２次元ランダムウォーク |
/* s0620.c
* ２次元ランダムウォーク
* (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;
}
|
拡散過程 (文献（１０）参照)
たくさんの粒子を一ケ所に集めて、ブラウン運動をさせると、粒子はどんどん広がって、ついには、一様にひろがってしまうでしょう。このような過程を拡散過程といいます。本節では、この拡散過程をシミュレーションします。
ところで、単位時間に１度ずつ、衝撃を受けるブラウン運動で、 t=0 で原点にいた粒子が、 t=n にはどこにいる確率が高いでしょうか。図６．３．１のように、一回の衝撃で動く距離を１とすると、ブラウン粒子が、 n 回目の衝撃をうけると、原点からの距離 dn は、つぎのように求められます。
 |
| ・・・・式６．３．１ |
となり、 n を大きくしていくと、平方根の中の２項目は n にくらべて無視できるようになります。したがって
 |
・・・・・ |
式６．３．２ |
となり、原点からの平均距離は衝撃の回収がふえると、その回数の平方根に比例して遠くなることがわかります。
たくさんの粒子を一ケ所に集めて、ブラウン運動をさせると、粒子はどんどん広がって、ついには一様に分布してしまいます。このことを、実際にシミュレーションで確かめてみましょう。
| 図６．３．１ 粒子の位置 |
 |
プログラムのポイント
プログラムs0630.cにおいては、２つのウインドウを用います。ウインドウ０では、粒子が運動する平面、ウインドウ１は、原点からの粒子の平均距離を表示することとします。
時刻 0 で、２０００個の粒子は原点にあるとし、クロックを１進めるたびに、前節と同じように、 0〜2π の乱数を発生させ、ランダムウォークをさせます。２０００個の粒子に対しこのことを実行し、原点からの粒子の平均距離
d を計算したあと、時刻を進めていきます。
実行結果
このプログラムを実行してみましょう。図６．３．２のように、はじめ原点にあった粒子は、あたかも生き物のように動きまわり、全体はしだいに拡散してい
く様子が表示されます。また、たしかに原点からの平均距離は、ウインドウ１に示すように、時刻の平方根に比例していることも示されます。
| 図６．３．２ 拡散のシミュレーション結果(途中の図） |
 |
時刻=10
|
 |
時刻=100
|
 |
時刻＝500
|
| s0630.c 拡散のシミュレーション（長時間ＣＰＵ占有に注意） |
/* 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 という枠をはみだしてしまいます。そこで、今回は枠まで粒子が到着したら、そこに壁があり、その壁ではねかえるようにします。
また、前節では、粒子が１クロックで動く長さを一定としていましたが、中心部では速く、周辺部では遅く動くようにしてみましょう。どんなことがおこるでしょうか 。
プログラムのポイント
プログラムs0640.cにおいては、１０００個の粒子をブラウン運動させます。ただし、初期状態として一様に分布させることとし、また粒子の動き方は、１クロックで進む長さを
l とすれば、
 |
・・・・・ |
式６．４．１ |
とします（ MAX は壁までの距離、 xi, yi は粒子 i の位置）。すなわち、中心部で最も速く、
 |
・・・・・ |
式６．４．２ |
でうごき、中心部から最も遠い隅の
 |
・・・・・ |
式６．４．３ |
の時、停止状態となるものとします。
また、次のクロックにおける位置を計算した結果、その値が MAX をこえた場合、図６．４．１のように、壁のところで、完全反射するものとします。すなわち
if (x[i] >= MAX) { x[i] = MAX - (x[i] - MAX) ; }
のように MAX からはみだした分、ひき返すように計算します。 y 座標も同様です。
| 図６．４．１ 粒子の動き方 |
 |
実行結果
このプログラムを実行してみましょう。部屋のほこりは、中心部では風が吹いたり、人の動きがあるため、速く動きまわり、次第に風の弱い部屋の隅に集まってきます。このような様子が図６．４．２のように、シミュレーションによって確かめることができます。
図６．４．２ ほこりはすみにたまるシミュレーション結果
開始時には均一に分布していたホコリは
時間がたつと下のように隅にたまっていく |
開始時
|
時刻=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 ６章１．解答 ２次元ランダムウォーク （歩幅がランダムに変わる場合） |
/* a0610.c ６章１．解答
* ２次元ランダムウォーク （歩幅がランダムに変わる場合）
* (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 が小さい場合のベルヌーイ試行をおこない、事象が起ったときに再度確率１／２で、右に行くか左に行くかきめるとします。このときの、軌跡を図示するプログラムを作りなさい。（歩幅一定、あるく時間が不定のランダムウォーク）
（解答）
| a0620.c ６章２．解答 １次元ランダムウォーク （歩幅一定，時間不定） |
/* a0620.c ６章２．解答
* １次元ランダムウォーク （歩幅一定，時間不定）
* (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;
}
|
|
 |
