[Matlab]如何用蒙地卡羅方法(Monte Carlo method)求出圓周率(π)?

最近在整理碩博時期的資料時，無意間看到博二的書報課程繳交的書面報告。該學期的書報課程有一個主題是電腦模擬在醫學物理上的應用，由長庚大學醫放系趙自強教授主講。對於他的演講內容，我已經非常模糊了。但印象最深刻的是，他明確提到蒙地卡羅方法(Monte Carlo method)常用於各領域的模擬計算，其中用此方法可以計算出圓周率π，也引起我寫一個小程式來驗證一下。

其原理的說明如上圖所示，假設有一個圓半徑為1，所以四分之一圓面積就為PI，而包括此四分之一圓的正方形面積就為1。如果隨意地將石頭丟入正方形中，則這些石頭有些會落於四分之一圓內，假設總共進去的石頭有n顆，落在圓內的石頭有c點，則依比例來算，就會得到圓周率PI的計算式了。至於如何判斷所產生的點落於圓內，只要令亂數產生X與Y兩個數值，如果X^2+Y^2等於1就是落在圓內。利用上述原理，要寫成程式碼就較為簡單了，利用matlab簡單驗證一下，以其內建亂數函數rand()當成任意丟出石頭的顆數，inside表示落入四分之一圓內的石頭顆數，最後圓周率PI可以計算得出。當然，以概率來說，隨著丟出的石頭顆數越多，其計算出來的圓周率也會越來越精確。真有趣~