シミュレーションをしよう(1):拡散

拡散のシミュレーション(1)

コップの水にインクを一滴落とすと、かきまぜなくてもインクはじわじわと水の中にひろがっていき、最終的は水とインクは均一にまじる(塩のようにイオンになって溶けるのではなく、細かい粒が水の中に浮んだ状態。コロイドという)。このような現象を拡散という。ミクロに見ると、ランダムに揺れ動く水分子がインクの粒と衝突することによってインクの粒がランダムに動く”ブラウン運動”が拡散の原因である。ここでは、粒子のランダムな運動をシミュレーションし、拡散現象を調べてみよう。日本語では「酔っ払いの歩き方」にたとえて酔歩と呼ばれる

一個

一個の粒子が一次元の線上をブラウン運動する様子をシミュレーションすることを考えてみよう。 たとえば、-1,0,1の三つの値をランダムに発生させて、出た数だけ粒子の位置を動かすことにする(-1なら左へ1歩、1なら右へ1歩動かし、0ならその場に留まる)。つまり、粒子の位置は整数値で表わされることになる。 これを実現するには次のようなプログラムが考えられる
#include<stdio.h>
#include "rnd.h"
#define MAXSTEP 100000
main()
{
  int i,ichi,step,randinit;
  int r[MAXSTEP], rseed[RND_P];

  scanf("%d%d",&step,&randinit);

  ichi = 0;
  ranset(randinit, rseed);

  rnd(r,step,rseed);
  for(i=0; i<step; ++i){
    ichi += (r[i]%3-1);
    printf("%d\n",ichi);
  }
}
ただし、乱数を発生するransetとrndは以前使ったものと同じである。 これは非常に簡単なプログラムで、 端末からステップ数と乱数の初期値を入力すると ichiの値をstep回計算して、そのたびに出力する。 初期値やステップ数を変えてみて、結果をgnuplotで図にしてみるべきだろう。

 見慣れない構文は下から2文目の
    ichi += (r[i]%3-1);
だろう。%記号は整数同士の割り算をした”余り”を与える演算子である。 r[i]%3は、r[i]の値を3で割った余りであり、値は0,1,2のいずれかに なる。つまり、この文は「i番目の乱数を3で割った余りから1を引いた数を ichiに足せ」という意味。なぜこれでよいか、考えてみること

N個

我々が現実に目にする拡散現象は膨大な数の微粒子が拡散するものである。 そこで、一個ではなく、N個の粒子を同時にブラウン運動させるようにプログラムを 変更してみよう。たとえば
#include<stdio.h>
#include "rnd.h"
#define N 20
main()
{
  int ichi[N];
  int step,randinit;
  int i,j;
  int r[N], rseed[RND_P];

  scanf("%d%d",&step,&randinit);

  for(i=0; i<N; ++i)ichi[i] = 0;
  ranset(randinit, rseed);

  for(j=1; j<step; ++j){
   rnd(r,N,rseed);
   for(i=0; i<N; ++i){
    ichi[i] += (r[i]%3-1);
   }
   for(i=0; i<N; ++i) printf(" %d",ichi[i]);
   printf("\n");
  }
}
これでいいだろう。 ただしこのプログラムは(Nかけるstep)個の数値を出力するので、 Nやstepを大きくすると出力量が多くなる(最近のコンピュータの 能力や容量を考えると、たいした問題ではないのだが)

分布

たとえば、粒子が100000個あって、1000ステップ後の粒子の分布を知りたいとしよう。すると、上のプログラムのままでは10万個の数が出力されてしまう。一方、1000ステップ後には最も遠くまで動く粒子でも左右に1000ずつしか進まないのだから、"どの場所に何個の粒子がいるか"を出力したほうが都合がよい。 下の例は、10000個の粒子でstep回計算後の粒子の分布を出力するものである。
#include
#include "rnd.h"
#define MAXSTEP 10000
#define N 10000
main()
{
  int ichi[N];
  int step,randinit;
  int i,j,k;
  int r[N], rseed[RND_P];
  int basho[2*MAXSTEP+1];

  scanf("%d%d",&step,&randinit);

  for(i=0; i<N; ++i)ichi[i] = 0;
  ranset(randinit, rseed);

  for(j=1; j<step; ++j){
   rnd(r,N,rseed);
   for(i=0; i<N; ++i){
    ichi[i] += (r[i]%3-1);
   }
  }

  for(k=0; k<2*MAXSTEP+1; ++k) basho[k] = 0;
  for(i=0; i<N; ++i) ++basho[ichi[i]+MAXSTEP];
  for(k=0; k<2*MAXSTEP+1; ++k) printf("%d %d\n",k-MAXSTEP,basho[k]);
}


プログラム中では下から2文目の
++basho[ichi[i]+MAXSTEP];
が一番わかりにくいだろう。何をしているか理解できるだろうか。

全体にもう少し工夫してみたいところではある
次へ