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

[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.5点：1/12日まで)

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

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

1：上下にのみ1づつ移動する場合を考える．上と下に移動する確率がそれぞれ1/2として計算を行う．n=100の場合について1000ステップ計算し，アニメを観察せよ．

2：移動した結果がマイナスとなる場合は0に戻るとしてn=100の場合について1000ステップ計算し，原点からの平均距離が時刻の何乗に比例して変化するかグラフを書くことで求めよ．

3：10ステップに一回，乱数で移動した後に必ず下に1移動するとする．原点からの平均距離が時刻の何乗に比例して変化するかグラフを書くことで求めよ．

4：1000ステップの計算が終了したのち，x=0, 1, 2, 3,...の各場所にいる球の数を求め，横軸x　縦軸　球の数のグラフを書く．球の数はxに対してどのような依存性を持っているか考察せよ．n=10000とする．アニメを作るためのwrite文は削除した上で実行すること．

日程表へ戻る << 　

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