モンテカルロ法で円周率を計算して、グラフにプロットする：無償ソフトで技術計算しよう【入門編】（2）（2/2 ページ）

» 2013年12月04日 10時00分 公開
[伊藤孝宏，MONOist]
前のページへ 1|2 　　　　　　

内部定数

　ヘルプに記載されている内部定数の一覧を下記に示します。これらは、所与の値として利用できます。例えば、「pi*5*5」とすると、半径5の円の面積を計算できます。また、「a=2+3*i」とすると、「a」は実数部「2」、虚数部「3」の複素数となります。上述のようにこれらは変数としても使うことが可能で、しかも警告なしで上書きされます。

  • Base Constants
  • E Euler Constant (Base of Natural Logarithm)：自然対数の底、e=2.718281828459045
  • EPS Double Precision Floating Point Relative Machine Precision Epsilon：倍精度での計算機誤差未満の数、eps　約2.2204e-16
  • FALSE Logical False：false=0
  • FEPS Single Precision Floating Point Relative Machine Precision Epsilon：単精度での計算機誤差未満の数、feps　約5.9604e-8
  • I-J Square Root of Negative One：i,j　虚数単位
  • INF Infinity Constant：inf　無限大と見なせる数、
  • PI Constant Pi：円周率、pi=3.14159265358979
  • TEPS Type-based Epsilon Calculation
  • TRUE Logical TRUE：true=1


　特に、「i」と「j」はCやJavaなどを使っている方にとって、forループ内のカウンターとして使いたくなりますが、別な変数名とするようにしましょう。これらに加えて、計算結果を示す変数ansと、数値ではないという意味の「NaN」（Not a Number）も変数として使うと上書きされるので注意が必要です。

例題：モンテカルロ法で円周率を計算する

　例題として、上述の円周率にちなんで、モンテカルロ法で円周率を計算してみましょう。モンテカルロ法とは、乱数を用いてシミュレーションを行う手法で、円周率を計算するには、x座標、y座標に0〜1までの乱数を設定した多数の点を生成して、原点からの距離が1以下の範囲にある点の数を数え、これと全体の点の数との比の4倍を求めます。

　FreeMatでは、下記のように2行でモンテカルロ法による円周率の計算が可能です。

--> x=rand(1000,1)+i*rand(1000,1);
--> length(find(abs(x)<1))/length(x)*4
ans = 3.1400

　詳しく説明すると、「rand(1000,1)」で「1000行1列の0以上1未満の乱数」を発生させます。変数xとして、「rand(1000,1)+i*rand(1000,1)」により「複素平面上に1000個の点を発生」させます（「複素平面」とは、実数部がX座標、虚数部がy座標の平面）。

　次に、複素数の大きさを計算します。複素数の大きさを計算するには、共役複素数との積、つまり、虚数部の符号が反転した複素数と掛け合わせた数を計算し、その平方根を求めます。例えば、「1+i」の大きさは「sqrt((1+i)*(1-i))」としますが、「abs(1+i)」でも計算できます。

　さて、大きさが1未満の点の数を選び出すには、「find(abs(x)<1)」とします（大きさが1以下ではなく1未満なのは、発生させた乱数が1未満であるためです）。こうすると、「x」という配列の中で大きさが「1」未満の要素の位置を返します。lengthは配列の行と列の内、大きい方のサイズを返します。例えば、3行2列の配列のlengthは「3」です。結局「length(find(abs(x)<1))」は大きさが1未満の点の数を示し、これを全体の点の数「length(x)」で割ると、点の数が十分多ければ、半径1の4分の1円と正方形との面積の比、すなわち「π/4」と等しいと見なせます。最後に4を掛けると円周率が求まります。

　円周率を求めるだけでは芸がないので、発生した点を図に表示させてみた結果が図2です。

図2　モンテカルロ法による円周率計算

　図2を表示させるのが下記の「ex01.m」というプログラムで、最初の3行で円周率を計算しています。つづく4〜6行で、複素平面状に大きさが1未満では青い点、1以上では赤い点となるようにプロットしています。8行目で半径1の4分の1円を描き、9、10行目でグラフの形を整えています。これらは入門編の後のグラフィックス編で説明します。

clear;N=1000;
x=rand(N,1)+i*rand(N,1);
length(find(abs(x)<1))/N*4
xin=x(find(abs(x)<1));
xout=x(find(abs(x)>=1));
plot(real(xin),imag(xin),'b.',real(xout),imag(xout),'r.');
hold('on');
plot(cos(deg2rad(0:90)),sin(deg2rad(0:90)),'k:');
axis('square');
axis([0,1,0,1]);
ex01.m

＞＞「ex01.m」ダウンロード

　上の「ex01.m」ファイルは、カレントフォルダあるいはパスを設定したフォルダに保存した後、「-->ex01」と入力すると、自動で計算してグラフを表示してくれます。エディタから実行させる場合は、左から2番目のアイコンをクリック、メニューの「File→Open File」、「Ctrl＋O」のいずれかでex01.mファイルを開き、エディタの左から10番目の右向き三角のアイコンをクリックするか、メニューの「Debug→Execute Current Buffer」とするか、F5キーを押すのいずれでも実行してくれます。

◇


　次回は演算子、関数について説明します。

参考文献

  • 「MATLABハンドブック」小林一行著、秀和システム刊
  • 「はじめてのFreeMat」赤間世紀著、工学社刊

無償ソフトで技術計算しよう

FreeMat

無償の工学計算ソフトでも、かなり高度な計算ができる！


筆者紹介

伊藤孝宏（いとう・たかひろ）

1960年生。小型モーターメーカーのエンジニア。博士（工学）。専門は流体工学、音・振動工学。現在は、LabVIEWを使って、音不良の計測・診断ソフト、特性自動検査装置などの開発を行っている。




前のページへ 1|2 　　　　　　

Copyright © ITmedia, Inc. All Rights Reserved.