2022年11月20日日曜日

FORTRNで一様乱数生成

euler.md

FORTRANで乱数生成

一様乱数を生成してみた。


方法

標準の "RANDOM_NUMBER" を使う。

REAL(8) の a を 0~1 の乱数にしたいなら、CALL RANDOM_NUMBER( a ) とする。


やってみる

PROGRAM RND REAL(8) :: a DO i = 1, 950000 ! 95万コの乱数生成 CALL RANDOM_NUMBER( a ) WRITE(100,*) a END DO END PROGRAM RND

コンパイル・実行

ifort -O2 rnd.f90 ./a.out

fort.100 というファイルができる。

これを見てもよくわからないので、ヒストグラムを作る。

PROGRAM HISTGRAM IMPLICIT NONE INTEGER, PARAMETER :: NUMBER=950000 ! データ個数 INTEGER, PARAMETER :: DIV=100 ! x軸の分割数 REAL(8) :: a(NUMBER), d(0:DIV) ! データと幅 INTEGER :: N(DIV), i, j ! 個数 d(0)=0.d0 ; N(:)=0 ! 初期化 ! データ読み込み OPEN( 100, FILE="fort.100", STATUS="OLD" ) DO i = 1, NUMBER READ(100,*) a(i) END DO CLOSE(100) ! ヒストグラムの幅を決める DO j = 1, DIV d(j) = d(j-1) + 1.d0/dble(DIV) END DO ! ヒストグラムの作成 DO i = 1, NUMBER DO j = 0, DIV-1 IF( a(i).GE.d(j) .AND. a(i).LT.d(j+1) ) N(j)=N(j)+1 END DO END DO ! 結果出力 >> fort.101 DO j = 0, DIV-1 WRITE(101,'(ES16.2, i6)') (d(j)+d(j+1))/2.d0, N(j) END DO END PROGRAM HISTGRAM

結果

gnuplotでヒストグラムは次のようにプロットする。

plot "fort.101" u 1:2 w boxes

結果は以下の通り。

例えば、0-1以外で乱数が欲しいときは。

CALL RANDOM_NUMBER( a ) a = a*10.d0 ! 0~10まで a = a -5.d0 ! -5~5まで

0 件のコメント:

コメントを投稿