LAMIS Sommerkonferanse 2025

Buffons Needle - Bergning av $\pi$

Manuelt fra målinger

Gitt:

• Et uendelig stort art
• Arket har parallelle linjer med avstand $d$
• Pinner med lengde $d$ slippes ned på arket

Utfall:

• Antall pinner som slippes: $N$
• Antall pinner som krysser en linje: $K$

Det kan vises at forholdet mellom antall pinner som slippes delt på antall pinner som krysser en linje, blir $\pi$ når du slipper veldig mange pinner.

$$ \lim_{N \to \infty} \frac{N}{K} = \pi $$

N = 22
K = 7

print("Beregner forholdet mellom N og K")
print(f"N / K = {N / K}")

La oss samle inn data

import numpy as np

N = [22, 25, 30, 28, 27, 40, 15]
K = [7, 8, 9, 8, 15, 9, 3]

print("Like mange oppføringer i listene?")
print(len(N)==len(K))

N = np.array(N)
K = np.array(K)
estimat = N / K
for i in range(len(N)):
    print(f"\nGruppe {i+1}: {estimat[i]:.4f}")

N_total = sum(N)
K_total = sum(K)

print("\nBeregner forholdet mellom N og K")
print(f"N_total / K_total = {N_total / K_total}")

Simulering

Her simuleres kast med $N$ nåler. Det er også mulig å justere avstanden $D$ mellom linjene og nålenes lengde $L$.

I simuleringen velges en tilfeldig $x$- og $y$-retning, samt en tilfeldig vinkel $\theta$.

Forklaring av de essensielle linjene:

crosses = np.any(((y1 < lines_y) & (lines_y < y2)) | ((y2 < lines_y) & (lines_y < y1)))

Som bygger på:

lines_y = np.arange(0, 10 + D, D)

og beregning av endepunkter basert på midtpunktet $(x_0, y_0)$ og vinkelen $\theta$.

Gjennomgang

1. lines_y: Dette er en liste med y-koordinatene til alle de horisontale linjene.

Hvis D = 2 (avstand mellom linjene), blir: lines_y = [0, 2, 4, 6, 8, 10]

2. y1 og y2: Dette er y-koordinatene til de to endene på nålen. De regnes ut fra midtpunktet (x0, y0) og vinkelen theta.

3. (y1 < lines_y) & (lines_y < y2) Her sjekker vi for hver linje (i lines) om den ligger mellom endepunktene y1 og y2.

Hvis y1 = 1.5 og y2 = 3.0, vil linjen på y=2 være mellom dem (True), men linjen på y=0 eller y=4 vil ikke (False).

4. (y2 < lines_y) & (lines_y < y1) Dette gjør akkurat det samme som over, men for tilfeller der y2 er mindre enn y1. (Det er bare for å ta høyde for at nålen kan helle nedover eller oppover.)

5. |: er som or, men er en elementvis operator. Den sammenligner hvert element i array 1 med det tilsvarende elementet i array 2, og gir tilbake et nytt array med resultatene. Den ser på to sammenligninger og gir True dersom en av testene gir True.

6. np.any(...) Dette sjekker om minst én av linjene ble krysset (om vi fikk en True i listen).

• Hvis nålen krysser minst én linje → crosses = True
• Hvis ikke → crosses = False.

Eksempel:

Avstand mellom linjer: D = 2 → lines_y = [0, 2, 4, 6, 8, 10] Nålen har endepunkter:

• y1 = 1.5
• y2 = 3.0

Sjekk: Linje på y=2 → ligger mellom endepunktene 1.5 og 3.0 → True Alle andre linjer → ligger ikke mellom → False

Resultat: np.any(...) finner minst én True → crosses = True (nålen krysser).

Beregning av $\pi$

Formelen som gir estimatet for $\pi$ er gitt ved:

$$ \pi_{\text{est}} = \frac{2 \cdot L \cdot N}{K \cdot D} $$

der:

• $L$ = lengden på nålen
• $D$ = avstand mellom linjene
• $N$ = totalt antall kast
• $K$ = antall kryssinger (nålen treffer en linje)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# --- Parametre -------------------------------------------------------------
D = 1.0        # Avstand mellom linjene
L = D / 2      # Lengde på nålene (L ≤ D)
N = 100         # Antall nåler som kastes

# --- Forbered figuren ------------------------------------------------------
x_min, x_max = -1, 11
y_min, y_max = -1, 7

lines_y = np.arange(0, 10 + D, D)          # y-koordinatene til linjene

plt.figure(figsize=(10, 6))
for y in lines_y:
    plt.plot([x_min, x_max], [y, y], 'k-', alpha=0.5, lw=1)

# --- Kast nålene -----------------------------------------------------------
cross_count = 0
segments = []                               # (x1, y1, x2, y2, crosses)

for _ in range(N):
    # Tilfeldig midtpunkt og vinkel
    x0 = np.random.uniform(0, 10)
    y0 = np.random.uniform(0, 6)
    theta = np.random.uniform(0, np.pi)

    # Endepunkter
    dx = (L / 2) * np.cos(theta)
    dy = (L / 2) * np.sin(theta)
    x1, y1 = x0 - dx, y0 - dy
    x2, y2 = x0 + dx, y0 + dy

    # Sjekk om nåla krysser en av linjene
    crosses = np.any(((y1 < lines_y) & (lines_y < y2)) |
                     ((y2 < lines_y) & (lines_y < y1)))

    segments.append((x1, y1, x2, y2, crosses))
    cross_count += int(crosses)

# --- π-estimat -------------------------------------------------------------
pi_est = (2 * L * N) / (cross_count * D) if cross_count else np.nan

# --- Plot nålene -----------------------------------------------------------
for x1, y1, x2, y2, crosses in segments:
    plt.plot([x1, x2], [y1, y2],
             color='red' if crosses else 'blue',
             lw=2, alpha=0.75)

# --- Pynt ------------------------------------------------------------------
plt.title(f"Buffons Needle  –  π ≈ {pi_est:.4f}  (Antall={N} / Krysser={cross_count})")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
#plt.xlabel("x")
#plt.ylabel("y")

plt.legend(handles=[
    mpatches.Patch(color='red',  label='Krysser'),
    mpatches.Patch(color='blue', label='Krysser ikke')],
    loc='upper right')

plt.tight_layout()
plt.show()

Simulering med mange nåler. Uten plotting

Enkel simulering

import numpy as np

D = 1.0
L = D / 2         # L ≤ D
N = 10_000

# 1) avstand fra senter til nærmeste linje. y0 er nå en array med N oppføringer.
y0 = np.random.uniform(0, D/2, N)

# 2) vinkel 0–π/2 (alltid positiv sin)
theta = np.random.uniform(0, np.pi/2, N)

# 3) Vi trenger kun å se på y-verdien, så ingen x-verdi er nødvendig
dy = (L/2) * np.sin(theta)

# 4) krysser linja?
crosses = y0 <= dy
k = crosses.sum()

pi_estimat = 2 * L * N / (D * k)
print(f"π ≈ {pi_estimat:.6f}  (Kryss: {k} av {N})")

Enkel simulering, ulike størrelsesordner

import numpy as np
from time import time

# Velg størrelsesorden
n = 8

t0 = time()
for i in range(1,n+1):
    
    D = 1
    L = D / 2
    N = 10**i
    
    y0 = np.random.uniform(0, D, N)
    theta = np.random.uniform(0, np.pi/2, N)
    dy = (L / 2) * np.sin(theta)
    y1 = y0 - dy
    y2 = y0 + dy
    
    kryss = (y1 < 0) | (y2 > D)
    antall_kryss = np.sum(kryss)
    pi_estimat = 2 * L * N / (D * antall_kryss)
    
    print(f"Estimert pi = {pi_estimat:.9f} med {antall_kryss:>6} kryssende nåler av {N} kast.")
print(f"\nFerdig! Det tok {(time()-t0):.2g} sekunder.")

Estimert pi = 5.000000000 	 med      2 kryssende nåler av 10 kast.
Estimert pi = 4.166666667 	 med     24 kryssende nåler av 100 kast.
Estimert pi = 3.125000000 	 med    320 kryssende nåler av 1000 kast.
Estimert pi = 3.126954346 	 med   3198 kryssende nåler av 10000 kast.
Estimert pi = 3.152883312 	 med  31717 kryssende nåler av 100000 kast.
Estimert pi = 3.150439801 	 med 317416 kryssende nåler av 1000000 kast.
Estimert pi = 3.140601921 	 med 3184103 kryssende nåler av 10000000 kast.
Estimert pi = 3.141489186 	 med 31832037 kryssende nåler av 100000000 kast.

Ferdig! Det tok 1.8 sekunder.