使用python实现布丰投针法

对于π的值，直到1946年的时候，人类才能将π的值精确计算到小数点后2037位，而现在的超级计算机的能力可以精确的计算到小数点后几十亿位，然而在计算机发明之前，还是使用这里的布丰投针法来计算π值，是最实用的方法。

使用代码来模拟这个过程，首先是程序设计思路的基本路程：


import  math
import  numpy as np
from  numpy import random
import  matplotlib.pyplot as plt
import  matplotlib.lines as mlines

然后根据代码设计的流程图实现代码：

UPDATE_FREQ = 500  #每模拟投针500次后就将结果绘制出来
BOUND = 10  #图形绘制区域的长和宽
BORDER = 0.05 * 10  #绘制区域边缘宽度
NEEDLES=  10000  #投针的次数
NEEDLE_LENGTH = 1  #针的长度
FLOORBOARD_WIDTH = 2  #两条平行线之间的距离
FLOORBOARD_COLOR = 'black'  #两条平行线的颜色
NEEDLE_INTERSECTING_COLOR = 'red'  #铁针与线相交时的颜色
NEEDLE_NON_INTERSECTING_COLOR = 'green'  #铁针与直线不相交时的颜色
class Needle:
    def  __init__(self, x = None, y = None, theta = None, length = NEEDLE_LENGTH):
        '''
        x, y 作为铁针的中点，由于图像绘制区域的宽度为BOUND，因此我们可以把中点坐标设置在宽度内的任何地方,
        theta是铁针与水平方向的夹角，在分析时我们将该值限制在0和Pi/2之间，考虑这个区间是因为我们确定哪个平行线
        与铁针跟接近的情况下。如果不考虑这个前提条件，那么夹角的值就可以限定在0和pi之间。我们模拟时不事先考虑那条线与
        针的距离更近，因此我们对夹角区间采用0和pi之间
        '''
        if x is None:
            x = random.uniform(0, BOUND)
        if y is None:
            y = random.uniform(0, BOUND)
        if theta is None:
            theta = random.uniform(0, math.pi)
        self.center = np.array([x, y])  #设置铁针中心点坐标
        self.comp = np.array([length/2 * math.cos(theta), length/2 * math.sin(theta)])  #根据铁针与水平方向的夹角，计算中心在水平方向和竖直方向的距离
        self.endPoints = np.array([np.add(self.center, -1 * self.comp),
                                  np.add(self.center, self.comp)])  #根据中心与水平方向和竖直方向的距离运算铁针两头的坐标
    def  intersectsY(self, y):
        return  self.endPoints[0][1] < y and self.endPoints[1][1] > y  #y是平行线在竖直方向上的坐标
class Buffon_Sim:  #启动模拟进程
    def  __init__(self):
        self.floorboards = []  #存储平行线在竖直方向上的y坐标
        self.boards = int ((BOUND / FLOORBOARD_WIDTH) + 1)  #计算平行线在绘制区域内的数量
        self.needles = []  #存储模拟的铁针
        self.intersections = 0  #记录铁针与平行线相交的数量
        window = "Buffon"
        title = "Buffon Needle Simulation"
        desc = (str(NEEDLES) + " needles of length " + str(NEEDLE_LENGTH) + 
               " uniformly distributed over a " + str(BOUND) + " by " + str(BOUND) + 
               " area" + " with floorboards of width " + str(FLOORBOARD_WIDTH))  #描述当前模拟情况
        self.fig = plt.figure(figsize = (8,8))
        self.fig.canvas.set_window_title(window)
        self.fig.suptitle(title, size = 16, ha = 'center')
        self.buffon = plt.subplot()  #将模拟投针绘制出来
        self.buffon.set_title(desc, style = 'italic', size = 9, pad = 5)
        self.results_text = self.fig.text(0, 0, self.updateResults(), size = 10)  #将投针情况绘制成图像
        self.buffon.set_xlim(0 - BORDER, BOUND + BORDER)
        self.buffon.set_ylim(0 - BORDER, BOUND + BORDER)
        plt.gca().set_aspect('equal')  
    def  plotFloorboards(self):
        for j in range(self.boards):
            self.floorboards.append(0 + j * FLOORBOARD_WIDTH)  #绘制平行线
            self.buffon.hlines(y = self.floorboards[j], xmin = 0,
                              xmax = BOUND, color = FLOORBOARD_COLOR, linestyle = '--', linewidth = 2.0)
    def  tossNeedle(self):  #模拟投针过程
        needle = Needle()
        self.needles.append(needle)
        p1 = [needle.endPoints[0][0], needle.endPoints[1][0]]  #获取铁针两个端点的x坐标
        p2 = [needle.endPoints[0][1], needle.endPoints[1][1]]  #获取铁针两个端点的y坐标
        for k in range(self.boards):  #检测铁针是否与平行线相交
            if needle.intersectsY(self.floorboards[k]):
                self.intersections += 1
                self.buffon.plot(p1, p2, color = NEEDLE_INTERSECTING_COLOR, linewidth = 0.5)  #将相交的铁针用红色线条表示
                return 
        self.buffon.plot(p1, p2, color = NEEDLE_NON_INTERSECTING_COLOR, linewidth = 0.5)  #不相交的铁针用绿色线条表示
       
        
    def  plotNeedles(self):
        
        for i in range(NEEDLES):
            self.tossNeedle()
            self.results_text.set_text(self.updateResults(i + 1))
            if (i + 1) % UPDATE_FREQ == 0:  #连续模拟指定次数投针后把结果绘制出来
                plt.pause(0.0001)     
    def  updateResults(self, needlesToTossed = 0):
        if self.intersections == 0:
            sim_pi = 0
        else :
            sim_pi = (2 * NEEDLE_LENGTH * needlesToTossed) / (FLOORBOARD_WIDTH * self.intersections)  #根据公式（4）计算Pi值
        error = abs(((math.pi - sim_pi) / math.pi) * 100)  #计算模拟结果与真实结果的误差
        s = ("Intersections: " + str(self.intersections) + \
               "\nTotal Needles: " + str(needlesToTossed) + \
               "\nApproximation of pi: " + str(sim_pi) + "\nError: " + str(error) + "%")
        return s
    def  plot(self):
        legend_lines = [mlines.Line2D([], [], color = FLOORBOARD_COLOR, linestyle = '--', lw = 2),
                       mlines.Line2D([], [], color = NEEDLE_INTERSECTING_COLOR, lw = 1),
                       mlines.Line2D([], [], color = NEEDLE_NON_INTERSECTING_COLOR, lw = 1)]
        self.buffon.legend(legend_lines, ['floorboard', 'intersecting needle', 'non-intersecting needle'],
                          loc = 1, framealpha = 0.9)
        self.plotFloorboards()
         
        self.plotNeedles()
        plt.show()
      
        
bsim = Buffon_Sim()
bsim.plot()

使用的公式以及代码运行的结果：