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

[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）にどう近づくか？

[2]ランダムウォーク

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

このようなプロセスは「ランダムウォーク」（酔歩）と呼ばれている．分子の拡散，熱の伝導など地球科学においても非常に重要なプロセスを表現している．

parameter文の利用

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

	integer x(100),y(100)

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

例：

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

練習問題3

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

注意点

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

例：

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

練習問題4

アニメーションを観察する．以下のテキストを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上で

gnuplot> n=1
gnuplot> load 'anim.txt'

とすることでアニメーションを観察すること．

練習問題5

球は原点からじょじょに遠ざかる．つまり，「拡散」する．格子上ではなくて連続空間を考えると，単位面積あたりに存在する球の数nは，次の拡散方程式に従って変化する．

$\begin{displaymath} {\partial n\over \partial t}=D\left({\partial^2n\over \partial x^2} + {\partial^2n\over \partial y^2}\right) \end{displaymath}$

(2)

拡散方程式は偏微分方程式であるが，微分を無理矢理割り算で置き換えると

$\begin{displaymath} {n\over t}=4D{n\over L^2} \end{displaymath}$

(4)

と書くことができる．ここでLは球が存在している範囲の距離である．これをtについて解くと

$\begin{displaymath} t={L^2\over 4D} \end{displaymath}$

(5)

となり，拡散した距離Lは時間の平方根に比例することがわかる．つまり，等速で運動するよりもはるかに遅い．

原点からの距離の平均値を各時刻で求め，時刻とともにファイルに出力する．n=10000とせよ．なお，アニメを観察するための write(1,*) x(i), y(i)は削除すること．巨大なファイルが作成されてしまう．
gnuplotで横軸時刻，縦軸距離の平均値で両対数グラフを作成する．x**0.5 と一緒にプロットすることで，時間の1/2乗に比例して平均距離が増加することを確かめる．

今日の宿題 (1/11日まで)

練習問題５のプログラム，および時刻ー平均距離のデータファイルを送付すること．

[6] 追加の宿題 (1/11日まで)

それぞれの球の原点からの距離(ここではdisとする）を整数型変数に代入する(たとえばndis=dis）．するとndisは0,1,2,...の値をとることになる．

ndisが0である球の数をnr(0), 1である球の数をnr(1)...といった具合に，それぞれの整数の距離にある球の数をカウントするための配列nr(0:1000)を用意する．

球の数を100000とし，1000ステップ計算する．各ステップにおいてnr(i)を算出する．各球のndisを出し，nr(ndis)=nr(ndis)+1で球の数を増やす．各ステップの開始前に，nr=0で数をリセットすること．

課題1: 時間とともに，nr(0)がどう変化するのか調べる．横軸時間，縦軸nr(0)のグラフを作り，時間の何乗に比例してnr(0)が変化しているのか求めよ．

課題2: 10, 100, 1000ステップ終了後の nr(i)をグラフにする．横軸i, 縦軸nr(i)のグラフを3つ作り，gnuplotで重ねてプロットすることで(たとえば plot '10.dat','100.dat','1000.dat')，球の存在している範囲が広がっていることを確認すること．

課題１の答えとプログラムを送付すること．

日程表へ戻る << 　

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