1. 程式人生 > >各種求圓周率π的演算法(蒙特卡洛法的Java實現)

各種求圓周率π的演算法(蒙特卡洛法的Java實現)

什麼是演算法?簡單地說,演算法就是有窮規則構成的用於解決某一類問題的運算序列或執行步驟。在《演算法之美:隱匿在資料結構背後的原理》第1章中我們講到要解決一個問題可能會有不同的方法,當時所舉的例子就是求圓周率π的近似值。對於這個問題你能想到多少種演算法呢?

探祕演算法世界,求索資料結構之道;

彙集經典問題,暢享程式設計技法之趣;

點撥求職熱點,敲開業界名企之門。

方法(1):首先根據微積分中的泰勒公式,可以使用下面這個公式來求π的近似值


方法(2):其次可以還可以利用下面這個馬青公式來求解π的近似值:


馬青公式由英國天文學家約翰·馬青(John Machin,1686 –1751)於1706年發現,他利用這個公式計算到了100位的圓周率。此外,上式中的反三角函式可以領用下面這個泰來展開式進行計算:


方法(3):最後,也是我們準備來動手實踐一下的就是利用蒙特卡洛法。這個方法背後的理論基礎是概率論中的伯努利大數定律。眾所周知,圓的面積公式為 πr2,當 r = 1 時,S = π。如果在一個邊長為1的正方形中有一個內切圓,那麼該圓的面積就是π/4。利用這一原理,可以假設有大量的隨機點等概率均勻地落入此正方形中。若正方形的面積為S',內切圓的面積為S,如果落入正方形中的點數為 N',落入內切圓的面積為N,這時便有 S'/S = N'/N,所以可據此求得π = N/N'

下面就編寫一個Java程式來實現蒙特卡洛法求π的近似值:

  1. import java.io.BufferedReader;  
  2. import java.io.InputStreamReader;  
  3. publicclass MonteCarloPi {  
  4.     publicstaticvoid main(String[] args) throws Exception{  
  5.         BufferedReader br = new BufferedReader(new InputStreamReader(System.in));  
  6.         System.out.print("How many darts shoud I toss at the board?\n");  
  7.         String s = br.readLine();  
  8.         int numberOfDarts = Integer.parseInt(s.trim());  
  9.         double radius = 1.0;  
  10.         Dartboard d = new Dartboard(radius);  
  11.         for(int i=1; i<=numberOfDarts; i++){  
  12.             Toss t = Toss.getRandom(radius);  
  13.             d.strike(t);  
  14.         }  
  15.         double fractionIn = d.getFractionIn();  
  16.         double pi = 4.0 * fractionIn;  
  17.         System.out.println("Pi is approximately "+pi);  
  18.     }  
  19. }  
  20. class Dartboard{  
  21.     privatedouble radius;  
  22.     privateint insideCircle, outsideCircle;  
  23.     public Dartboard(double radius){  
  24.         this.radius = radius;  
  25.         insideCircle = 0;  
  26.         outsideCircle = 0;  
  27.     }  
  28.     publicvoid strike(Toss toss){  
  29.         double x = toss.getX();  
  30.         double y = toss.getY();  
  31.         if(Math.sqrt(x*x + y*y) < radius)  
  32.             insideCircle++;  
  33.         else
  34.             outsideCircle++;  
  35.     }  
  36.     publicdouble getFractionIn() {  
  37.         double total = (double) (insideCircle + outsideCircle);  
  38.         return (double) insideCircle/total;  
  39.     }  
  40. }  
  41. class Toss{  
  42.     privatedouble x,y;  
  43.     public Toss(double x, double y){  
  44.         this.x = x;  
  45.         this.y = y;  
  46.     }  
  47.     publicdouble getX(){return x;}  
  48.     publicdouble getY(){return y;}  
  49.     publicstatic Toss getRandom(double radius){  
  50.         double x,y;  
  51.         double size = Math.random();  
  52.         double sign = Math.random();  
  53.         size = size * radius;  
  54.         if(sign > 0.5)  
  55.             x = size;  
  56.         else
  57.             x = -size;  
  58.         size = Math.random();  
  59.         sign = Math.random();  
  60.         size = size * radius;  
  61.         if(sign > 0.5)  
  62.             y = size;  
  63.         else
  64.             y = -size;  
  65.         returnnew Toss(x,y);  
  66.     }  
  67. }  

執行一下上述程式碼,然後我們來驗證一下其效果如何...
  1. How many darts shoud I toss at the board?  
  2. 1000000  
  3. Pi is approximately 3.142028  

可見所得到的結果較好地近似於圓周率 π 。

(本文完)