布丰投针：用随机实验估算π值

创作时间:

作者:

@小白创作中心

布丰投针：用随机实验估算π值

引用

CSDN

https://blog.csdn.net/Miles_zhg/article/details/145919964

当我们观察这个世界时，会发现一些令人惊叹的现象，其中布丰投针问题就是一个典型例子。通过随机投掷针的方法，竟然能够估算出π的值，这一发现令人着迷。下面将详细介绍布丰投针问题的原理及其计算方法。

问题描述

布丰投针问题是18世纪由Georges-Louis Leclerc, Comte de Buffon首次提出的。假设我们有一个由平行木条组成的地板，每条木条的宽度都是a，我们将一个长度为L（L < a）的针随意扔到地板上，那么针与木条相交的概率是多少？

数学推导

设x为针的中点到最近一条线的距离，0 ≤ x ≤ a/2；φ为针与线的角度。参考图（a）

全集Ω = { (φ, x) | 0 ≤ φ ≤ π, 0 ≤ x ≤ a/2 }。

距离x和角度φ相同的情况视为同一种情况。全集为夹角π和a/2围成的坐标。参考图（b）

相交的条件为：针的一半作为斜边，x/sinφ ≤ L/2。

不相交的条件为：x/sinφ > L/2。

事件G = { (φ, x) | 0 ≤ φ ≤ π, 0 ≤ x ≤ L/2 sinφ }。

其中，L为针的长度，a为两条平行线的距离。

定积分计算概率

P = ∫(0 to π) (L/2) sinφ dφ / (π * a/2) = (L/2) [-cosπ + cos0] / (πa/2) = L/(πa/2) = 2L/(πa)

从公式可以看出，针越长（L越大），概率越大；平行线距离越大（a越大），概率越小。

蒙特卡罗方法

可以通过概率求出π。例如，当a为L的两倍时，有π = 1/P。

Python实现

import numpy as np

def buffon(n, l=0.8, a=1):
    k = 0
    theta = np.random.uniform(0, np.pi, n)
    x = np.random.uniform(0, a / 2, n)
    
    for i in range(n):
        if x[i] <= (l * np.sin(theta[i])) / 2:
            k += 1
    
    if k == 0:
        return np.inf
    else:
        return (2 * l * n) / (a * k)

print(buffon(10**8))

输出结果：

3.14187290067662

Plotly仿真

import numpy as np
import plotly.graph_objects as go
import plotly.io as pio

def generate_needles(num_needles, length=1, line_spacing=2):
    # 生成随机的针数据
    x_centers = np.random.uniform(0, 5, num_needles)
    y_centers = np.random.uniform(0, 3 * line_spacing, num_needles)
    thetas = np.random.uniform(0, np.pi, num_needles)
    # 计算到最近平行线的距离
    closest_lines = np.round(y_centers / line_spacing) * line_spacing
    distances = np.abs(y_centers - closest_lines)
    # 判断是否相交
    intersects = distances <= (length / 2) * np.sin(thetas)
    # 计算端点坐标
    half_length = length / 2
    x0 = x_centers - half_length * np.cos(thetas)
    y0 = y_centers - half_length * np.sin(thetas)
    x1 = x_centers + half_length * np.cos(thetas)
    y1 = y_centers + half_length * np.sin(thetas)
    return [
        {
            "x0": x0[i],
            "y0": y0[i],
            "x1": x1[i],
            "y1": y1[i],
            "intersects": intersects[i],
        }
        for i in range(num_needles)
    ]

# 创建基础图形
fig = go.Figure()
# 添加平行线
line_spacing = 2
y_lines = np.arange(0, 3 * line_spacing + line_spacing, line_spacing)
for y in y_lines:
    fig.add_hline(y=y, line=dict(color="green", width=2))
# 预生成所有针数的数据并创建帧
frames = []
steps = []
for n in range(1, 1001, 50):
    needles = generate_needles(n)
    red_x, red_y = [], []
    blue_x, blue_y = [], []
    for needle in needles:
        if needle["intersects"]:
            red_x += [needle["x0"], needle["x1"], None]
            red_y += [needle["y0"], needle["y1"], None]
        else:
            blue_x += [needle["x0"], needle["x1"], None]
            blue_y += [needle["y0"], needle["y1"], None]
    # 创建帧
    frame = go.Frame(
        data=[
            go.Scatter(x=red_x, y=red_y, mode="lines", line=dict(color="red")),
            go.Scatter(x=blue_x, y=blue_y, mode="lines", line=dict(color="blue")),
        ],
        name=str(n),
    )
    frames.append(frame)
    # 创建滑块步骤
    steps.append(
        {
            "args": [[frame.name], {"frame": {"duration": 0}, "mode": "immediate"}],
            "label": str(n),
            "method": "animate",
        }
    )
# 添加初始数据（空数据）
fig.add_scatter(x=[], y=[], mode="lines", line=dict(color="red"), name="相交")
fig.add_scatter(x=[], y=[], mode="lines", line=dict(color="blue"), name="未相交")
# 配置动画和滑块
fig.frames = frames
fig.update_layout(
    title="布丰投针模拟实验 (拖动滑块后点击▶️播放)",
    xaxis=dict(range=[0, 5], showgrid=False),
    yaxis=dict(range=[0, 6], showgrid=False),
    updatemenus=[
        {
            "type": "buttons",
            "buttons": [
                {
                    "args": [None, {"frame": {"duration": 0}}],
                    "label": "▶️",
                    "method": "animate",
                }
            ],
        }
    ],
    sliders=[{"active": 0, "currentvalue": {"prefix": "针数: "}, "steps": steps}],
    width=800,
    height=600,
)
# 生成HTML文件，嵌入plotly.js
pio.write_html(
    fig, file="buffon_simulation.html", auto_play=False, include_plotlyjs=True
)  # 改为True确保本地嵌入
print("已生成 buffon_simulation.html，请用浏览器打开")