Keith Vetter 2003-02-18 - a Monte Carlo simulation of Buffon's Needle experiment.
In 1777, Georges Louis Leclerc, Comte de Buffon [
1] asked and solved the following problem: what is the probability that a needle of length L when dropped onto a wooden floor with evenly spaced cracks, distance D apart, will land on a crack?
The answer, P = 2L/pD, leads to an amusing way to estimate the value of Pi.
One amusing
gotcha in doing the simulation that trips up a lot of people is this: to simulate dropping the needle, one random value you need is the angle the needle is dropped at. If you're not careful, you might code this as
set angle [expr 2*$PI*rand()]. You're using Pi to estimate Pi! (If you use degrees you're probably just letting the language library do the conversion to radians for you.) A cute way around this problem is to exploit the fact that Pi is transcendental and just use the consecutive values 1,2,3....
KBK Hmm, another cute way to draw a random angle is:
- Choose a random number x, -1 <= x <= 1.
- Choose another random number y, -1 <= y <= 1
- If x**2 + y**2 > 1, go back to the top.
- (x,y) is now a random point inside the unit circle. Return atan2( y,x )
This method doesn't use Pi to estimate Pi, but does choose a random angle, unlike the method of simply using integers.
KPV One side benefit of simply using integers is that it can be much faster. If you exploit the fact that
sin(n+1) = sin n A· cos(1) + cos n A· sin(1), and
cos(n+1) = cos n A· cos(1) - sin n A· sin(1), you can compute the sine value in just 4 floating point multiplies and 2 additions. This will eventually get round off errors but not for a long while.
##+###################################################
#
# buffon.tcl -- Simulation of Buffon's Needle Problem
# by Keith Vetter - Feb 9, 2003
#
set S(unit) 40 ;# Length of needle
set S(num) 500
set S(tm) 0
set S(draw) 0
set S(total) 0
set S(cross) 0
set S(cosn) 1 ;# We'll incrementally...
set S(sinn) 0 ;# compute sin and cos
set S(cos1) [expr {cos(1)}]
set S(sin1) [expr {sin(1)}]
proc DoDisplay {} {
wm title . "Buffon's Needle"
pack [frame .ctrl -relief ridge -bd 2 -padx 5 -pady 5] \
-side right -fill both -ipady 5
canvas .c -relief raised -bd 2 -height 500 -width 500 -bg \#ffe4c4
.c xview moveto 0 ; .c yview moveto 0
pack .c -side top -fill both -expand 1
DoCtrlFrame
for {set y $::S(tm)} {$y < 4000} {incr y $::S(unit)} {
.c create line 0 $y 4000 $y -tag grid
}
Clear
update
trace variable ::S(draw) w Tracer
set ::S(draw) 0
}
proc DoCtrlFrame {} {
scale .s -from 1 -to 5000 -orient h -variable S(num) -label Needles \
-relief ridge -highlightthickness 0
button .start -text "Start" -command Start -bd 4
button .clear -text "Clear" -command Clear -bd 4
button .stop -text "Stop" -command {set S(draw) 0} -bd 4
.start configure -font "[font actual [.start cget -font]] -weight bold"
.clear configure -font [.start cget -font]
.stop configure -font [.start cget -font]
frame .pi -bd 2 -relief raised -pady 5
label .lbl -text "PI Estimation" -bd 2 -relief groove
label .ltot -text "Total"
label .etot -textvariable S(total) -width 6 -bd 2 -relief sunken -bg white \
-anchor w
label .lcross -text "Cross"
label .ecross -textvariable S(cross) -width 6 -bd 2 -relief sunken -bg white \
-anchor w
label .lpi -text "PI"
label .epi -textvariable S(pi) -width 6 -bd 2 -relief sunken -bg white \
-anchor w
button .about -text About -command About
grid .s - -in .ctrl -sticky ew -row 0
grid .start - -in .ctrl -sticky ew -row 2
grid .clear - -in .ctrl -sticky ew
grid .stop - -in .ctrl -sticky ew -pady 10
grid .pi - -in .ctrl -sticky ew -row 11
grid .lbl - -in .pi -sticky ew -row 0
grid .ltot .etot -in .pi -sticky ew -row 2
grid .lcross .ecross -in .pi -sticky ew
grid .lpi .epi -in .pi -sticky ew
grid rowconfigure .ctrl 1 -minsize 10
grid rowconfigure .ctrl 10 -minsize 50
grid rowconfigure .ctrl 50 -weight 1
grid columnconfigure .ctrl 1 -weight 1
grid rowconfigure .pi 1 -minsize 10
grid .about - -in .ctrl -row 100 -sticky ew
}
proc Clear {} {
.c delete needle
set ::S(total) 0
set ::S(cross) 0
set ::S(pi) ""
}
proc NextSin {} {
global S
set s [expr {$S(sinn) * $S(cos1) + $S(cosn) * $S(sin1)}]
set c [expr {$S(cosn) * $S(cos1) - $S(sinn) * $S(sin1)}]
set S(sinn) $s
set S(cosn) $c
return [list [expr {$S(unit) * $s}] [expr {$S(unit) * $c}]]
}
proc Start {} {
set ::S(draw) 1
DropNeedle $::S(num)
}
proc About {} {
set msg "Buffon's Needle\nby Keith Vetter, February 2003\n\n"
append msg "A Monte Carlo method to estimate the value of pi invented\n"
append msg "Comte de Buffon in 1777.\n\n"
append msg "If you drop a 3 inch needle onto a wooden floor with 3 inch\n"
append msg "slats, the probability that the needle will cross a crack\n"
append msg "is 2/pi. Unfortunately, it is very poor way to get an\n"
append msg "accurate value of pi.\n"
tk_messageBox -title "About Buffon's Needle" -message $msg
}
proc Tracer {var1 var2 op} {
if {$::S(draw) == 0} { ;# Stopping drawing
.stop config -state disabled
.start config -state normal
.clear config -state normal
} else { ;# Starting to draw
.stop config -state normal
.start config -state disabled
.clear config -state disabled
}
}
proc DropNeedle {{cnt 1}} {
global S
for {set i 0} {$i < $cnt} {incr i} {
if {$S(draw) == 0} break
.c itemconfig needle -width 1
set w [winfo width .c]
set h [winfo height .c]
# Drop needle between lines visible on the screen
set x [expr {$w * rand()}]
set maxy [expr {int(($h - $S(tm))/$S(unit)) * $S(unit)}]
set y [expr {$S(tm) + $maxy * rand()}]
# Get other end point
foreach {s c} [NextSin] break
set x1 [expr {$x + $c}]
set y1 [expr {$y + $s}]
set yy [expr {($y - $S(tm))}]
set yy [expr {$yy - $S(unit) * int($yy / $S(unit))}]
set yy [expr {$yy + $s}]
set color black
incr S(total)
if {$yy <= 0 || $yy >= $S(unit)} { ;# Did it cross
set color red
incr S(cross)
}
set S(pi) [expr {$S(cross) / .2 / $S(total)}]
.c create line $x $y $x1 $y1 -fill $color -tag needle -width 5
if {$i < 50 || ($S(total) % 50) == 0} update
}
set S(draw) 0
}
DoDisplay
uniquename 2013jul29
This code could use an image to show what it produces:
(Thanks to 'gnome-screenshot', 'mtpaint', and ImageMagick 'convert' on Linux for, respectively, capturing the screen to a PNG file, cropping the image, and converting the resulting PNG file to a JPEG file that was less than half the size of the PNG. Thanks to FOSS developers everywhere --- including Linux kernel and Gnu developers. I used the 'mv' command and the ImageMagick 'identify' command in a shell script to automagically rename the cropped file to its current pixel dimensions.)
This GUI comes up with no needles shown --- only horizontal lines. I clicked on the 'Start' button to generate the default of 500 needles. Note that the needles that lie across a horizontal line are shown in red. The needles that lie between the lines show as black.
