Function Repository Resource:

GaltonBoardModel

Source Notebook

Model ball trajectories in a 2D Galton board using Hertzian force laws between the ball and the pegs

Contributed by: Michael Trott

ResourceFunction["GaltonBoardModel"][brd,inits, output]

models the dynamics in a Galton board brd with ball initial conditions inits and returns the model result output as the graphic or data.

ResourceFunction["GaltonBoardModel"][brd,prop]

computes the property prop of the Galton board brd.

ResourceFunction["GaltonBoardModel"]["RandomGaltonBoard"]

returns a random Galton board.

Details and Options

A Galton board board is specified as an association with the following keys specifying the pegs, the rows of pegs and the ball parameters (all assumed to be dimensionless numbers):
"RowCount"number of rows of pegs
"PegRadius"radius of the pegs
"PegHardness"hardness of the pegs
"PegForceExponent"exponent of the force law
"PegRotation"rotation of the peg deformation from the radial direction
"VerticalAcceleration"vertical acceleration
"LinearFrictionCoefficient"friction force on the balls proportional to their speed
"QuadraticFrictionCoefficient"friction force on the balls proportional to their speed squared
"BallRadius"radius of the balls
WLOG, the position of the center of the single peg of the first row, is assumed to be {0,0}. The vertical spacing between rows is fixed to be 1.
For small deformations, the force law is assumed to be of elastic Hertzian form force~deformationα, where α is the exponent of the force law of the peg-ball interaction. (The sphere-sphere force law is α=3/2, and the cylinder-cylinder force law with aligned axis is α=1.) For large deformations, the force law smoothly approaches a maximum value.
The model numerically solves the Newtonian equations of motion for the balls until a ball reaches one unit below the last row of pegs.
Each ball is assumed to move through the Galton board individually; no ball-ball interactions are included.
The ball radius can be 0.
The initial conditions are lists, and each element can be of the following form:
xiassume y=1 and {vx,i,vy,i}={0,0}
{xi,yi}assume {vx,vy}={0,0}
{{xi,yi},{vx,i,vy,i}}fully specified initial conditions
The following output types output can be specified:
"PathGraphics"a graphic of the pegs and ball paths
"PathGraphics3D"a 3D graphic of the pegs and ball paths with the ball velocities as the z values
"Paths"list of paths as styled line primitives
"FinalPosition"a list of the positions of the balls one unit down from the last row of pegs
"PathArcLengths"a list of the arc lengths of the ball paths
"ArrivalTimes"a list of the arrival times of the ball paths
"InterpolatingFunctions"a list of pairs of {position,velocity} interpolating functions of the ball paths
{"PathArrayPlot", dim}an array plot of dimension dim×dim of the ball path
"TimeEvolution"an interactive demonstration allowing you to step through the ball-falling process
"Manipulate"an interactive demonstration allowing you to change board parameters and ball initial positions
For the "Manipulate" output type, at least one of the parameter values or initial conditions must be an interval.
The following options can be specified:
"TrajectoryStyle"{blue,white}color directive for the ball paths
"PegStyle"slightly opaque dark redpair of center-boundary colors of the pegs
The local color of a peg is a blend of the two "TrajectorStyle" colors, with the blending ratio according to the Hertzian force.
The following properties prop are supported for boards:
"TotalPegCount"the total count of pegs
"MinPegDistance"the minimum distance between two pegs
"PegPositions"a list of the center corrdinates of the pegs
"ForceLawSketch"sketch of the specified Hertzian force law

Examples

Basic Examples (5) 

A single ball path in a large Galton board of 60 rows of pegs:

In[1]:=
board = <|
   "RowCount" -> 60, "PegRadius" -> 0.25,
   "PegHardness" -> 100,
   "PegForceExponent" -> 1.25,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.4,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.1|>;
ResourceFunction["GaltonBoardModel"][board, {0.3}, "PathGraphics"]
Out[2]=

The ball paths for a typically sized Galton board using a realistic force law and finite-sized balls:

In[3]:=
board = <|
   "RowCount" -> 18, "PegRadius" -> 0.25, "PegHardness" -> 100,
   "PegForceExponent" -> 3/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.1|>;
ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 100], "PathGraphics"]
Out[3]=

Point-size balls, no friction lead and a small friction coefficient:

In[4]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> 0.3, "PegHardness" -> 100,
   "PegForceExponent" -> 3/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.5,
   "QuadraticFrictionCoefficient" -> 0,
   "BallRadius" -> 0.|>;
ResourceFunction["GaltonBoardModel"][board, Table[{x, 1}, {x, -1, 1, 2/48}], "PathGraphics"]
Out[4]=

Large, soft pegs allow particles to substantially penetrate and deform the pegs:

In[5]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> 0.4, "PegHardness" -> 3,
   "PegForceExponent" -> 2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.3,
   "QuadraticFrictionCoefficient" -> 0.2,
   "BallRadius" -> 0.05|>;
ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 100], "PathGraphics"]
Out[5]=

Sketch of the force law of α=1/2 and hardness 12 of a peg:

In[6]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> 0.4, "PegHardness" -> 12,
   "PegForceExponent" -> 1/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.3,
   "QuadraticFrictionCoefficient" -> 0.2,
   "BallRadius" -> 0.05|>;
ResourceFunction["GaltonBoardModel"][board, "ForceLawSketch"]
Out[6]=

The force law for α=4:

In[7]:=
ResourceFunction["GaltonBoardModel"][
 Append[board, "PegForceExponent" -> 4], "ForceLawSketch"]
Out[7]=

Scope (2) 

Use nonvanishing initial velocities for the balls (all originating from a single point):

In[8]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.32, "PegHardness" -> 400,
   "PegForceExponent" -> 1,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.|>;
ResourceFunction["GaltonBoardModel"][board, Table[{{0, 1}, 3 RandomReal[{-1, 1}, 2]}, 120], "PathGraphics"]
Out[8]=

A list of result types returns an association of results:

In[9]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> 0.5, "PegHardness" -> 20,
   "PegForceExponent" -> 3,
   "PegRotation" -> 1.2,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.4,
   "QuadraticFrictionCoefficient" -> 0.3,
   "BallRadius" -> 0.02|>;
ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 24],
                                      {"FinalPositions", "PathGraphics"}]
Out[10]=

PathSpeedGraphics3D (1) 

Show the speed of the balls as the vertical function value:

In[11]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.32, "PegHardness" -> 40,
   "PegForceExponent" -> 1,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.|>;
ResourceFunction["GaltonBoardModel"][board, Table[{{0, 1}, 3 RandomReal[{-1, 1}, 2]}, 60], "PathSpeedGraphics3D"]
Out[11]=

PathArrayPlot (1) 

Displaying many (tens of thousands or even hundreds of thousands) of paths would result in large graphics with many ball paths. Returning the result as an array plot is much more memory efficient:

In[12]:=
ResourceFunction["GaltonBoardModel"][<|
  "RowCount" -> 18, "PegRadius" -> 0.3, "PegHardness" -> 100,
  "PegForceExponent" -> 3/2,
  "PegRotation" -> 0,
  "VerticalAcceleration" -> 1,
  "LinearFrictionCoefficient" -> 0.8,
  "QuadraticFrictionCoefficient" -> 0,
  "BallRadius" -> 0.1|>, Table[{RandomReal[{-1, 1}], 1}, 100], {"PathArrayPlot", 1200}]
Out[12]=

TimeEvolution (1) 

Return an interactive demonstration/animation showing the time-dependent paths of the balls:

In[13]:=
ResourceFunction["GaltonBoardModel"][<|
  "RowCount" -> 12, "PegRadius" -> 0.3, "PegHardness" -> 20,
  "PegForceExponent" -> 0.8,
  "PegRotation" -> 0,
  "VerticalAcceleration" -> 1,
  "LinearFrictionCoefficient" -> 0.4,
  "QuadraticFrictionCoefficient" -> 0.2,
  "BallRadius" -> 0.1|>, Table[{x, 1}, {x, -1, 1, 2/8}], "TimeEvolution"]
Out[13]=

Manipulate (2) 

Interactively explore the ball paths as a function of some of the board parameters:

In[14]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> 0.4, "PegHardness" -> Interval[{1, 10.}],
   "PegForceExponent" -> 2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> Interval[{0, 4}],
   "QuadraticFrictionCoefficient" -> Interval[{0, 4}],
   "BallRadius" -> Interval[{0, 0.3}]|>;
ResourceFunction["GaltonBoardModel"][board, Table[{x, 1}, {x, -1, 1, 2/5}], "Manipulate"]
Out[14]=

Also include the initial conditions in the interactive controls:

In[15]:=
board = <|
   "RowCount" -> 12, "PegRadius" -> Interval[{0, 0.6}], "PegHardness" -> Interval[{1, 10.}],
   "PegForceExponent" -> 2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.2,
   "QuadraticFrictionCoefficient" -> 0.2,
   "BallRadius" -> Interval[{0, 0.3}]|>;
ResourceFunction[
 "GaltonBoardModel"][board, {Interval[{0, 2}], Interval[{0, 2}]}, "Manipulate"]
Out[15]=

ArrivalTimes (4) 

The ball paths in a relatively bouncy set of pegs:

In[16]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.25, "PegHardness" -> 80,
   "PegForceExponent" -> 2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.5,
   "QuadraticFrictionCoefficient" -> 0.1,
   "BallRadius" -> 0.1|>;
ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 20], "PathGraphics"]
Out[17]=

Compute the arrival (at the imaginary row 25) for 1,000 balls:

In[18]:=
 tAs = ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 1000], "ArrivalTimes"];

Plot the histogram of the distribution of the arrival times:

In[19]:=
hg = Histogram[tAs, {1}, "PDF", PlotRange -> All]
Out[19]=

Fit with a typical first-hit distribution for Brownian motion with drift:

In[20]:=
bins = {#1, #2/1000} & @@@ Tally[Round[tAs]];
NonlinearModelFit[bins, a Exp[-(b - c t)^2/(2 t)]/Sqrt[2 Pi t^3], {{a, 50}, {b, 42}, {c, 1}}, t] // Normal
Out[21]=
In[22]:=
Show[{hg, Plot[%, {t, 20, 75}]}]
Out[22]=

InterpolatingFunctions (3) 

In modeling the result for a Galton board with large pegs, many balls never penetrate the interior regions between the pegs:

In[23]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.4, "PegHardness" -> 8,
   "PegForceExponent" -> 5/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0,
   "QuadraticFrictionCoefficient" -> 2,
   "BallRadius" -> 0.02|>;
SeedRandom[1234];
ResourceFunction["GaltonBoardModel"][board, Table[{{0, 1}, 10 Normalize[RandomReal[{-1, 1}, 2]]}, 10], "PathGraphics"]
Out[24]=

Obtain the interpolating functions for the positions and velocities for 10 paths:

In[25]:=
SeedRandom[1234];
ifs = ResourceFunction["GaltonBoardModel"][board, Table[{{0, 1}, 10 Normalize[RandomReal[{-1, 1}, 2]]}, 10], "InterpolatingFunctions"];
Short[ifs, 2]
Out[26]=

Show the speed versus the vertical coordinate for all paths:

In[27]:=
yV[t_Real, {ifX_, ifV_}] := {Last[ifX[t]], Norm[ifV[t]]}
ParametricPlot[Evaluate[Style[yV[t, #], RandomColor[]] & /@ ifs],
 {t, 0, 20}, AspectRatio -> 1/2, Frame -> True, Axes -> False, FrameLabel -> {"y", "|\!\(\*OverscriptBox[\(V\), \(\[RightVector]\)]\)|"},
 PlotRange -> {All, {0, 1}}]
Out[27]=

Options (2) 

Specify yellow pegs and blue ball paths:

In[28]:=
board = <|
   "RowCount" -> 16, "PegRadius" -> 0.25, "PegHardness" -> 100,
   "PegRotation" -> 0,
   "PegForceExponent" -> 1,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.5,
   "QuadraticFrictionCoefficient" -> 0.3,
   "BallRadius" -> 0.|>;
ResourceFunction["GaltonBoardModel"][board, Table[{x, 1}, {x, -1, 1, 2/120}], "PathGraphics", "TrajectoryStyle" -> Directive[Thickness[0.001], Opacity[0.8], Lighter[Blue, 0.3]],
 "PegStyle" -> {Directive[Yellow, Opacity[0.66]], White}]
Out[28]=

To visualize a finite ball size, specify Thickness[Automatic] in the "TrajectoryStyle" option setting to obtain "thick" ball paths:

In[29]:=
board = <|
   "RowCount" -> 18, "PegRadius" -> 0.2, "PegHardness" -> 100,
   "PegRotation" -> 0,
   "PegForceExponent" -> 3/2,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.5,
   "QuadraticFrictionCoefficient" -> 0,
   "BallRadius" -> 0.2|>;
ResourceFunction["GaltonBoardModel"][board, Table[{x, 1}, {x, -1, 1, 2/36}],  "PathGraphics", "PegStyle" -> {Gray, White},
 "TrajectoryStyle" -> Directive[Thickness[Automatic], Opacity[0.2], Darker[Purple, 0.6]]]
Out[29]=

Applications (2) 

Show regression to the center for a board:

In[30]:=
board = <|
   "RowCount" -> 18, "PegRadius" -> 0.25, "PegHardness" -> 40,
   "PegForceExponent" -> 3/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.4,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.1|>;
{gr, fps} = ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 1000],
                                                             {"PathGraphics", "FinalPositions"}] // Normal;

The balls' paths:

In[31]:=
gr[[2]]
Out[31]=

The final positions of the balls:

In[32]:=
Short[fps[[2]], 3]
Out[32]=

Bin the data under the last row of pegs:

In[33]:=
binnedData = Round[First /@ fps[[2]], Sqrt[3]/2];

Show a histogram of the binned data:

In[34]:=
Histogram[binnedData]
Out[34]=

Fit the final positions to a normal distribution:

In[35]:=
estDist = EstimatedDistribution[binnedData, NormalDistribution[\[Mu], \[Sigma]]]
Out[35]=

Because of the binning, a direct test fails:

In[36]:=
DistributionFitTest[binnedData, estDist]
Out[36]=

Use bootstrapping and collect statistics:

In[37]:=
testStatisticSample = Block[{bootstrap, gd},
   Table[bootstrap = RandomVariate[estDist, Length[bins]];
                gd = Round[bootstrap, Sqrt[3]/2]; DistributionFitTest[gd, estDist, "TestStatistic"], {100}]];

Compute the p-value:

In[38]:=
NProbability[
 DistributionFitTest[binnedData, estDist, "TestStatistic"] > \[ScriptCapitalS], \[ScriptCapitalS] \[Distributed] testStatisticSample]
Out[38]=

Compute the autocorrelation function for the ball position.

Start the balls from deep within the board:

In[39]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.45, "PegHardness" -> 100,
   "PegForceExponent" -> 5/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.05|>;
ResourceFunction["GaltonBoardModel"][board, Table[{{0, -9}, 10 Normalize[RandomReal[{-1, 1}, 2]]}, 50], "PathGraphics"]
Out[39]=

Get the interpolating functions of the ball paths:

In[40]:=
ifs = ResourceFunction["GaltonBoardModel"][board, Table[{{0, -9}, 10 Normalize[RandomReal[{-1, 1}, 2]]}, 10], "InterpolatingFunctions"];

Define the position autocorrelation function by averaging along each path as well as over all paths:

In[41]:=
autoCorrelation[if_InterpolatingFunction, \[Tau]_, f_ : Identity] := Module[{T = if[[1, 1, 2]]}, Mean[Table[
    EuclideanDistance[f@if[t], f@if[t + \[Tau]]], {t, 0, T - \[Tau], \[Tau]}]]]
In[42]:=
positionAutoCorrelation[ifList_, \[Tau]_] := Mean[autoCorrelation[First[#], \[Tau]] & /@ ifList]

Plot c(τ) as a function of τ:

In[43]:=
pacData = Table[{\[Tau], positionAutoCorrelation[ifs, \[Tau]]}, {\[Tau], 0.1, 5, 0.1}];
ListLogLogPlot[pacData]
Out[43]=

Fit c(τ) to a power law c(τ)~τβ; the value β≲0.5 indicates subdiffusive movement:

In[44]:=
Fit[pacData, {1, \[Tau]}, \[Tau]]
Out[44]=

Properties and Relations (1) 

Very small vertical acceleration turns the Galton board into a classical Lorenz gas:

In[45]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.4, "PegHardness" -> 12,
   "PegForceExponent" -> 0.3,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 0.1,
   "LinearFrictionCoefficient" -> 0,
   "QuadraticFrictionCoefficient" -> 0.2,
   "BallRadius" -> 0.1|>;
ResourceFunction["GaltonBoardModel"][board, Table[{{0, -11}, 5 Normalize[RandomReal[{-1, 1}, 2]]}, 5], "PathGraphics", "PegStyle" -> {Gray, White},
  "TrajectoryStyle" :> Directive[Opacity[0.6], Thickness[0.002], RandomColor[]]] // Quiet
Out[45]=

Possible Issues (1) 

Negative friction coefficients (modeling autonomous balls) and rotating peg forces (modeling rotating pegs) can lead to unusual path shapes:

In[46]:=
board = <|
   "RowCount" -> 16, "PegRadius" -> 0.5, "PegHardness" -> 4,
   "PegForceExponent" -> 2,
   "PegRotation" -> 1.4,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> -2,
   "QuadraticFrictionCoefficient" -> 2,
   "BallRadius" -> 0.1|>;
ResourceFunction["GaltonBoardModel"][board, Table[{{0, -7}, 5 Normalize[RandomReal[{-1, 1}, 2]]}, 120], "PathGraphics", "PegStyle" -> {Gray, White},
 "TrajectoryStyle" :> Directive[Opacity[0.6], Thickness[0.002], RandomColor[]]]
Out[46]=

Neat Examples (3) 

Inject balls with slightly different initial conditions to see the trajectories diverge after multiple peg collisions:

In[47]:=
board = <|
   "RowCount" -> 24, "PegRadius" -> 0.5, "PegHardness" -> 12,
   "PegForceExponent" -> 1/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0.3,
   "QuadraticFrictionCoefficient" -> 0.2,
   "BallRadius" -> 0.05|>;
In[48]:=
ResourceFunction["GaltonBoardModel"][board, Table[{{2 + 0.4 RandomReal[{-1, 1}], 1},
                                                                    {-5, -4}}, {60}], "PathGraphics",
 "TrajectoryStyle" :> Directive[Opacity[0.6], Thickness[0.004], RandomColor[]],
 "PegStyle" -> {Directive[Gray, Opacity[0.5]], White}]
Out[48]=

Retrieve some "random" boards with preselected parameters:

In[49]:=
boards = Prepend[Evaluate[KeyDrop[#, "RowCount"]], "RowCount" -> 12] & /@ ResourceFunction["GaltonBoardModel"]["RandomGaltonBoard", 4]
Out[49]=
In[50]:=
TextGrid[Function[p, Prepend[#[p] & /@ boards, p]] /@ Keys[boards[[1]]], Dividers -> All, Alignment -> Left]
Out[50]=
In[51]:=
GraphicsRow[
 ResourceFunction["GaltonBoardModel"][#, Table[{x, 1}, {x, -1, 1, 2/24}], "PathGraphics"] & /@ boards]
Out[51]=

Connect ball trajectories and form an intricate polygon:

In[52]:=
board = <|
   "RowCount" -> 18, "PegRadius" -> 0.25, "PegHardness" -> 100,
   "PegForceExponent" -> 3/2,
   "PegRotation" -> 0,
   "VerticalAcceleration" -> 1,
   "LinearFrictionCoefficient" -> 0,
   "QuadraticFrictionCoefficient" -> 1,
   "BallRadius" -> 0.1|>;
In[53]:=
paths = ResourceFunction["GaltonBoardModel"][board, Table[{RandomReal[{-1, 1}], 1}, 48], "Paths"];
In[54]:=
Show[Graphics[
  Polygon[Join @@ MapIndexed[If[OddQ[#2[[1]]], #1, Reverse[#1]] &, First /@ Cases[paths, _Line, \[Infinity]]]]]]
Out[54]=

Version History

  • 2.0.0 – 11 January 2021
  • 1.0.0 – 06 November 2019

License Information