Galton Board

2025, Gabriel Lemieux

Resources

Wiki: https://en.wikipedia.org/wiki/Galton_board
​Matlab: https://github.com/minoue-xx/Galtonboard-simulator/tree/main
Real Life Demo: https://upload.wikimedia.org/wikipedia/commons/d/dc/Galton_box.webm
​

Display

display::usage="display[ball,pegs] Display the ball and the pegs.";​​display[pt_Association,pg_Association]:=Module[{board,histogram},​​board=Graphics[{​​(*DrawPegs*)​​{LightGray,EdgeForm[Gray],Disk[{#1,#2},pg[r]]}&@@@pg[pos],​​(*DrawBalls*)​​{EdgeForm[{Blue,Opacity[0.8]}],FaceForm[{Blue,Opacity[0.5]}],Disk[{#1,#2},pt[r]]}&@@@pt[pos]​​},​​PlotRange->{{-10,10},{-0,10}},​​PlotRangeClipping->True,​​ImageSize->Large,​​PlotLabel->StringForm["Galton Board (`` Balls)",Length[pt[pos]]+Length[pt[bins]]]​​];​​histogram=Histogram[pt[bins],{Range[-10,10]},"Probability",​​PlotRange->{{-10,10},Automatic},​​Axes->{False,True},AspectRatio->1/6,AxesOrigin->{0,0},AxesStyle->Directive[Gray,12],​​ImageSize->Large,ChartStyle->"Pastel"​​];​​Return@Rasterize[Column[{board,histogram}]]​​];​​(*Preventexecutionofilldefinedargument*)​​display[pt_,pg_]:=Failure["Bad Format, requires associations.",<|"pt"->Head[pt],"pg"->Head[pg]|>];​​

Simulation

(*Protectsymbolsusedaskeytopreventmistakes*)​​Protect[r,Δt,ϵ,g,pos,mom,bins];
(*Initialization*)​​ball=<|r->0.1(*Radiusoftheball*),n->100(*Numberofsimultaneousballs*)|>;​​peg=<|r->0.25(*Radiusofthepegs*)|>;​​sim=<|​​ Δt->0.05(*HigherValueneedtocorrectcollisionpositions*),​​ ϵ->0.4,(*1.0meanallenergybackintheball*)​​ g->{0,-9.81},(*Youcantilttheboard*)​​ n->1000(*TotalNumberofBalls*)​​|>;​​​​(*InitializePositionandMomentumMatrices*)​​randomBall:={RandomReal[{-0.05,0.05}],RandomReal[{0,1}]+10};​​ball[pos]=Table[randomBall,ball[n]];​​ball[mom]=Table[0,ball[n],2];​​ball[bins]={};​​​​(*Pegposition*)​​peg[pos]=Flatten[Table[{ii+0.5jj,jj},{ii,-15,10},{jj,-1,8}],1];​​​​(*DisplayLatestFrame(Quickerwithoutit)*)​​Dynamic[lastImage]​​​​(*MainLoop*)​​res={};​​While[Length[ball[pos]]>0,​​(*Step*)​​(*SchemeBackwardEuler*)​​(*NextStepMomentum*)​​ball[mom]+=Threaded[sim[Δt]sim[g]];​​​​(*Predict/Correctseemslowerthatdirectlycomputingthepos,butmoreaccuratesoallowhigherΔt*)​​(*NextPositionPrediction*)​​next=ball[pos]+sim[Δt]ball[mom];​​​​(*Distance:Donotparallelize,itsalreadyinthealgo.*)​​dist=DistanceMatrix[next,peg[pos]];​​​​(*Collisionsbetweencurrentandfuturestep*)​​col=Parallelize@MapIndexed[{v,i}|->{SelectFirst[v,(#<ball[r]+peg[r])&->"Index"],i[[1]]}/.({Missing[_],_}->Nothing),dist];​​​​(*CorrectMomentumtoincludecollisions*)​​If[Length[col]>1,​​{up,val}=Transpose[Parallelize@MapApply[{p,b}|->{b,sim[ϵ]*Norm[ball[mom][[b]]]Normalize[next[[b]]-peg[pos][[p]]]},col]];​​momt=ball[mom];​​momt[[up]]=val;​​ball[mom]=momt;​​];​​​​(*Updateposwithcorrectedmomentum*)​​ball[pos]+=sim[Δt]ball[mom];​​​​(*Filterthroughtheballs,remove/resetandregisterthoseexiting*)​​exited=Select[ball[pos],#[[2]]<0&->"Index"];​​ingame=Select[ball[pos],#[[2]]>0&->"Index"];​​If[Length[exited]>0,ball[bins]=Join[ball[bins],ball[pos][[exited,1]]]];​​​​(*Removethemfromtheposandmommatrices*)​​ball[pos]=ball[pos][[ingame]];​​ball[mom]=ball[mom][[ingame]];​​​​(*Insertnewball*)​​left=sim[n]-Length[ball[bins]]-Length[ball[pos]];​​toAdd=Min[left,ball[n]-Length[ball[pos]]];​​Do[AppendTo[ball[pos],randomBall];AppendTo[ball[mom],{0,0}],toAdd];​​​​(*Saveforanimation*)​​AppendTo[res,ball];​​​​(*Drawthelastimage*)​​lastImage=Rasterize@display[ball,peg];​​];​​(*Makethefinalvideo*)​​video=AnimationVideo[display[res[[i]],peg],{i,1,Length[res],1},DefaultDuration->sim[Δt]*Length[res]]