趣味の数学とプログラム。楽しい数学の世界へようこそ!

【Python】ビュフォンの針

ビュフォンの針(適当に投げた針が、平行線と交差する確率を求める問題)をPythonでシミュレーションしました。確率の求め方は、こちらの記事をご覧ください。

平行線の間隔を\(d\)、針の長さを\(l  (l\leq d)\)とし、針の始端Aのx座標、y座標は一様乱数で与えます。

針の終端Bを次のようにして求めます。

始点Aを中心とした単位円内にランダムに点を取り、Aからの距離をdx、dyとします。

三平方の定理から、r を求められるので、\(\sin\theta\)、\(\cos\theta\)を求めることができます。(\(\theta\)は、X軸と針の角度)

\(\sin\theta = \large{\frac{dy}{r}}\)、\(\cos\theta = \large{\frac{dx}{r}}\)

Aの位置を\((x,y)\)とすると、終端Bの位置は、\((x+l\large{\frac{dx}{r}} \normalsize{,y+l} \large{\frac{dy}{r}}\))

 

針が線と交差するのは

・AかBのy座標が、いずれかの平行線のy座標と等しい

・AとBのy座標をd(間隔)で割った商が異なる

のいずれかのときです。

 

平行線の間隔\(d\)を10、針の長さ\(l\)を5とし、200×200の平面上でN回(N本)試行するプログラムを掲載します。

 

'''
ビュホンの針で円周率計算(三角関数は使わない)
  平行線の間隔 10cm
  針の長さ   5cm
  200cm×200cm内で(N回 or N本)試行
'''

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

#  試行回数 or 本数の指定
N=1000  

#  線と交差する針 ゼロクリア
count_red=0

# 針の始端セット (一様乱数)
x=np.random.uniform(5.0,195.0,N) 
y=np.random.uniform(5.0,195.0,N)

#  グラフのサイズ指定
plt.figure(figsize=(10,10))   


# 全部の針の終端を計算 → 描画
for i in range(N):
  r=2                         # 単位円内にランダムに点を取る (0,0)は除く                         
  while r>1:
    dx=np.random.uniform(-1.0,1.0)
    dy=np.random.uniform(-1.0,1.0)
    if dx==0 and dy==0:       
      r=2
    else:
      r=math.sqrt(dx*dx+dy*dy) 

  x_end=x[i]+5*dx/r           #針の終端の座標
  y_end=y[i]+5*dy/r
  X=[x[i],x_end]              #針の座標
  Y=[y[i],y_end]
          
  if y[i]%10==0 or y_end%10==0: #端が線上 (赤)
    plt.plot(X,Y,color="r")
    count_red+=1
  else:
    if y[i]//10!=y_end//10:     #線と交差 (赤)
      plt.plot(X,Y,color="r")
      count_red+=1
    else:
      plt.plot(X,Y,color="c")   #線に触れない (シアン)
 
  
print("N = ",N)
print("線と交わる針:",count_red)
print("円周率の近似値:",N/count_red)


#平行線描画(黒)

X=[0,210]
for i in range(21):
  Y=[i*10,i*10]
  plt.plot(X,Y,color="k")



ax=plt.gca()
ax.xaxis.set_major_locator(ticker.MultipleLocator(10))  #X軸 目盛り指定
ax.yaxis.set_major_locator(ticker.MultipleLocator(10))  #Y軸 目盛り指定

plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です