数値解析法及び演習第十回

[1]乱数の発生法

数値実験を行うと乱数が必要になる場合がある．コンピュータで乱数を発生させる方法はいくつかあって，それについて学ぶ．なお，コンピュータで発生できるのは「真の」乱数ではなく疑似乱数である．

混合合同法

これは，

$\begin{displaymath} x_{n+1}=ax_n+b \quad {\rm mod}\,\mu \end{displaymath}$

(1)

によって順次乱数を発生させるものである．右辺で計算した値をμで割り，その余りを左辺に代入する．

練習問題1

3ビットの整数において乱数列を発生させ，その周期についてしらべてみよ． x₁=1とし，μ=4とせよ．a=3,b=1, a=2,b=1, a=1,b=3の場合についてしらべよ．

32ビット整数（我々が使っているinteger）では，混合合同法をもちいると

   implicit none
   integer i,j
   double precision ran
   
   do j = 1, 100
       i=843314861*i+453816693

      if(i.lt.0) i=(I+2147483647)+1
      ran=dfloat(i)*4.656612873D-10
     write(*,*) j,ran
    enddo
	
	end

をつかうことで[0,1)に規格化された乱数をつくることができる．

練習問題2

上のアルゴリズムを用いて乱数を発生させ，平均値，標準偏差を求めよ．乱数を発生させる回数を10, 100, 1000, 10000, 100000, 1000000と変化させてみよ．真の値（平均値=1/2, 標準偏差=1/(12)**0.5）にどう近づくか？

[2]ランダムウォーク

一辺が長さ1.D0の正方形からできている碁盤の目を考える．原点(0,0)にn個の球が置いてある．乱数を発生させ，上下左右に存在する隣の格子点にそれぞれ確率1/4で移動する．n回乱数を発生させ，n個全ての球を一回移動させると１ステップが終了する．

parameter文の利用

このプログラムを作成しようとすると，球の数nの大きさの配列を複数用意する必要がある．たとえばi番目の球のx座標，y座標．100個の場合なら

	double precision x(100),y(100)

などとと記述する．さらに，それぞれの球の原点からの距離もあったら dis(100)　などが追加される．同じ数字を何回も書くのは面倒であり，配列のサイズを変化させる場合に書き間違いが起こる．そこでparamtere文を用いる．parameter文は，ある文字列をその数字として読み替える，ということを行う．つまりプログラム中でその文字列が登場すると，parameter文で定義された数値にコンパイル時におきかわる．注意点は，parameter文で定義された変数は，プログラム中で別の値を代入はできない．ということである．

例：

      integer n
      parameter (n=1000)
      double precision x(n),y(n)

練習問題3

n個の球がランダムウォークするプログラムを作成せよ．始めはすべて原点に球は存在する．確率 1/4 = 0.25D0 で上下左右に1.D0だけ移動する．まずはn=10でプログラムを実行する．10000ステップ計算を実行すること．

注意点

n個のx座標，y座標をファイルに書き出したら，その次には空白行を入れておくこと．

例：

   do i=1,n
     write(1,*) x(i),y(i)
   enddo
     write(1,*)  ! 空白行

課題1

アニメーションを観察する．以下のテキストをanim.txtで保存する．

  set xrange [-100:100]
  set yrange [-100:100]

  if (n>1000) exit

  plot "ファイル名" every 1:1::n::n with points pt 7 ps 2

  pause 0.05

  n=n+1

  reread

gnuplot上で