3次元離散ランダムウォーク

暇を持てあましたhat-tunの遊び。

3次元離散ランダムウォークの到達地点距離のシミュレーションプログラムを書いてみた。
あとそのスケーリング。

3次元離散ランダムウォークってのは、
時間0の時に原点にいるヒトが、
各単位時間ごとにそれぞれ6分の1の確立で、x,y,zの軸の正負の方向に動き続けるというシンプルなもの。

それで、原点にわしゃわしゃいっぱいヒトがいた時に、何時間後にみんなだいたいどこにいるやろか。
って話ですね。

たとえば、実用例としては、
物理現象の、ブラウン運動などは3次元ランダムウォークですねー◎

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

#define N_SUM 100000 //試行回数
#define N_TOTAL 1000 //ステップ数上限

int main(void);
double frand(void);
void random_walk(void);


int main(){
  random_walk();
  return 0;
}


double frand(){
  return (double)rand()/RAND_MAX;
}


void random_walk(void){
  
  int i,j,k;
  int x,y,z;
  
  double r=0.0; //各試行での2乗和を累計
  double r_ave=0.0; //2乗和の平均
 
  FILE* fp;
  fp = fopen("test.dat","w+");
 
  for(i=1;i<=N_TOTAL;i=i*10){
    for(j=0;j<N_SUM;j++){
      
      x=0;
      y=0;
      z=0;
      srand((unsigned)time( NULL ));//乱数系列の更新      

      for(k=0;k<i;k++){
	double val = frand()*6;//0<=val<=6の乱数
	
	if(val>5.0){
	  x++;
	}else if(5.0>= val && val > 4.0){
	  x--;
	}else if(4.0>= val && val > 3.0){
	  y++;
	}else if(3.0>= val && val > 2.0){
	  y--;
	}else if(2.0>= val && val > 1.0){
	  z++;
	}else{
	  z--;
	}
	
      }
      
      r = r+ pow(abs(x),2)+pow(abs(y),2)+pow(abs(z),2);     
    }
    r_ave= (double)r/N_SUM;
    printf("%4dstep | R average = %f\n",i,r_ave);
    fprintf(fp,"%f %f\n",pow(i,0.5),pow(r_ave,0.5));
  }
  fclose(fp);
}

テキトーに書きました。上の例で言うと、N_SUMが人の数で、N_TOTALが時間ですね。
それにしてもステップ数のオーダーがO(N_SUM*N_TOTAL^2)というひどいプログラムですねw
まともに計算したら10^11となって、理論上1000秒もかかるみたいですw
ここでは、計算時間を1秒ぐらい(＝約1億ステップ)に抑えるため、
N_TOTALが1,10,100,1000の時の値を出してます。

それと乱数を2つ使ったほうが高精度という話を聞いたんですが、
よくわかりませぬ＞＜
教えて賢いヒト。

まぁそれはさておき、こっからが新しい。俺の中で。

gnuplotを使って、スケーリングを確かめまする。
上のプログラムで、
ステップ数、末端間距離の2乗平均=r_ave
の各1/2乗をファイルにアウトプットする。

test.datのoutput例

1.000000 1.000000
3.162278 3.710881
10.000000 16.904752
31.622777 28.647591

やりたいことは、この2つの変数が、大体比例してるやん！（スケーリングしてる）
ということを確認したいわけです。

そのためにgnuplotのfit関数というものを使うのですが、
こいつは、ある指定した（予想した）関数とデータから、その関数の係数を割り出してくれるものです。
さらにグラフにプロットすることもできます。
これで、データが予想した関数に近いかどうかを判定できるというスンポー。

論より証拠、やってみましょう。
とりあえずおもむろにgnuplotを起動する。
fittingする関数を定義する。　→　　机で計算して予想。
そしてデータファイルを読ませる。→　だいたいの係数がわかる
plotしてみる。　→　わっほい！
という流れ。

$ gnuplot
gnuplot> f(x)=a*x
gnuplot> fit f(x) 'test.dat' using 1:2 via a
gnuplot> plot f(x),'test.dat' using 1:2

ここでは、線形になるとの予想があるので、f(x)=a*xと定義。
そして、fit!
using 1:2というのはおまじない。いや、意味知らないだけですが。
via aで対象変数を設定しましょう。

せっかくなので、プロットしたグラフをファイルに。

gnuplot> set term png
gnuplot> set output "scaling.png"
gnuplot> plot f(x),'test.dat' using 1:2

と、ここではscaling.pngというファイルに出力してみます！
できたできた♪　　　どれどれどれ。。。

ん？
比例してねえwww

すいません。多分コードの精度の問題ですww
理論上は比例するはずです・・・w
ダメ過ぎるw
ああ、乱数まわりちょっと勉強せんとなー。orz

まぁこういう使い方があるってことだけでも、理解してくれたら幸い。