Keith Vetter 2003-02-18 - a Monte Carlo simulation of Buffon's Needle experiment.
In 1777, Georges Louis Leclerc, Comte de Buffon 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/�D, leads to an amusing way to estimate the value of pi.
On 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....
##+###################################################
#
# 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
}
DoDisplayCategory Graphics | Category Mathematics | Category Application | Category Whizzlet