ガウス分布の乱数を利用した速度の初期化

つづいて、2次元の面心格子上に配置された粒子に初期速度を与えるルーチンを作成します。今回は初期速度を与える手法として、正規分布（ガウス分布）を有する乱数を用います。その前に、C言語で一様（疑似）乱数を発生させる基本的なプログラムを下記に示します。

#include <stdio.h>
#include <stdlib.h>

main()
{
int i;
double r;
     srand(55);
     for(i = 0; i < 100; ++i){
          r = (double)rand()/RAND_MAX;
          printf("%lf\n", r);
     }
}

このプログラムでは、100個の乱数を発生・表示します。注意事項は、

• stdlib.hをインクルードすること
• srand()関数で乱数の系列を初期化
• rand()関数を呼び出すたびに0〜RAND_MAX(stdlib.hの中で定義されているのでそのまま使える)の間の疑似乱数整数が得られる
• 0〜1に規格化するためにRAND_MAXで割る場合、(double)rand()としてrand()の型を変換する必要があることに注意

では、続いて乱数の発生度数がある分布（例えばガウス分布）に従うような乱数が必要な場合はどうするかということですが、例えば、ヒットミス法と呼ばれる以下の方法を利用することで実現できます。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

main(){
  int i;
  double r, m, sigma;
  m = 0.5;
  sigma = 0.1;

  for(i = 0; i < 10000; ++i){
    while(1){
      r = (double)rand()/RAND_MAX;
      if(1.0/sqrt(2*M_PI*sigma)*exp(-(r-m)*(r-m)/2.0/sigma/sigma) >= (double)ran
d()/RAND_MAX)
        break;
    }
    printf("%lf\n",r);
  }
}

この場合、平均0.5、標準偏差0.1となる乱数10000個が発生されます。コンパイル後、例えば、0〜0.1の範囲にある乱数の数を調べるのには下記のコマンドが便利です。

./a.out | gawk '{if ($0 > 0.0 && $0 < 0.1) print $0}' | wc

なお、２つの一様乱数からガウス分布を有する２つの乱数を発生する方法にボックス・ミューラー法がある。こちらを調べて利用しても良い。

以上の練習をもとに、前節で作成した粒子に平均値が0.5となる標準偏差0.7となるガウス分布を有する速度($v_x, v_y$)を与えよ。これら速度は配列velocity[200][2]を用意しそれに格納すること。