Application Center - Maplesoft

# Stochastic Model of Filling a Box with Spheres

You can switch back to the summary page by clicking here.

SpheresInBox.mws

Stochastic Model of Filling a Box with Spheres

Lee R. Partin

lpartin@chartertn.net

February 11, 2002

The task is to determine how many balls will fit into a given size box during the random filling of the box.  This Maple document first calculates estimates to establish a lower and upper bound on the solution.  Then, a stochastic model is programmed to randomly place balls into the box.

The task was based upon the fizz ball example presented to me by Wayne Pafko.  ( http://www.pafko.com/wayne)

The box dimensions are set at 6-3/8" x 3" x 5-3/4".  The spheres have a radius of 3/4".

Lower and Upper Bounds on Solution

Box dimensions:

 > restart;

 > Length:=6+3/8: Width:=3: Height:=5+3/4:

Include sphere plotting routine,

 > with(plottools):

Lower Bound

Use Cubes Instead of Spheres for the Dimensions

Fill the box with cubes that have side lengths equal to 2 times the sphere radius.

 > Side:=2*Radius; CountX:=floor(Length/Side); CountY:=floor(Width/Side); CountZ:=floor(Height/Side); Count:=CountX*CountY*CountZ; # total cubes that fit into box LowerBound:=Count;

Upper Bound

Totally Fill the Volume

The number of spheres in the box can not be more than the amount of spheres that have a total volume equal to the box volume.

Tightly Packed Spheres

Spheres can pack tightly as shown in the following plot.  Note that only one sphere is shown for the middle layer in the plot.

In this packing method, there are two complete spheres within a box having sides with lengths of Layer3Height x 2*Radius x 2*Radius.  The volume taken up a one sphere is then:

If this packing method could fully use the box volume, the new upper bound estimate is:

 > UpperBound:=floor(BoxVolume/SphereSpace);

Source Code for the Stochastic Model of Filling the Box

The model randomly places equal size spheres into a box starting with the first layer followed by placing new spheres on top of existing ones.

Several geometic calculations are required for the model.  Maple is used to derive the formulas.  Then, a module called Box is programmed to bundle the calculations into a working package.

Geometric Model 1 (3 Spheres to Support New Sphere)

Given three spheres located in 3D, find the origin of a fourth sphere that rests on top of them.  Return FAIL for the z coordinate if it is not feasible.

Let:  r = radius of the spheres

(x1,y1,z1), (x2,y2,z2), (x3,y3,z3) = origin of the three given spheres

Then, the origin of the resting sphere is the point where the three given spheres would intersect if their diameters were doubled.  The algebra is solved by Maple as follows:

 > eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2: eqn3:=(x-x3)^2+(y-y3)^2+(z-z3)^2=(2*r)^2: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions.  They are converted to Maple functions as follows:

 > assign(solution): fx:=unapply(x,x1,y1,z1,x2,y2,z2,x3,y3,z3,r): fy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,y3,z3,r): fz:=unapply(z,x1,y1,z1,x2,y2,z2,x3,y3,z3,r):

The functions may return multiple values.  Also check for an infeasible solution:

 > fx2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;   try   res:=[allvalues(fx(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end: fy2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;   try   res:=[allvalues(fy(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end: fz2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;   try   res:=[allvalues(fz(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end:

Testing the results,

 > x1:=-0.5: y1:=-1: z1:=0:   # sphere 1 x2:=1.75: y2:=1.5: z2:=0:  # sphere 2 x3:=0: y3:=2.5: z3:=0:     # sphere 3                  r:=1:     # radius

 > x4:=fx2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r); y4:=fy2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r); z4:=max(fz2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r));

Viewing the Plot of Spheres

 > c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): c4:=sphere([x4,y4,z4], r): plots[display]([c1,c2,c3,c4], scaling=constrained,                style=patch, axes=boxed);

Geometric Model 2 (2 Spheres and Fixed x to Support New Sphere)

Given two spheres located in 3D and fixed x origin location for a third sphere, find the origin of the third sphere that rests on top of the two spheres.  Return FAIL for the z coordinate if it is not feasible.

Let:  r = radius of the spheres

(x1,y1,z1), (x2,y2,z2) = origins of the two given spheres
x3 = known x location of the resting sphere

Then, the origin of the resting sphere is the point where the two given spheres would intersect if their diameters were doubled.  The algebra is solved by Maple as follows:

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2: eqn3:=x=x3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions.  They are converted to Maple functions as follows:

 > assign(solution): gy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r): gz:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r):

The functions may return multiple values.  Also check for an infeasible solution:

 > gy2:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gy(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end: gz2:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gz(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end:

Testing the results,

 > x1:=1.25: y1:=1: z1:=0:   # sphere 1 x2:=2.25: y2:=3: z2:=0:   # sphere 2 x3:=1:                    # sphere 3 x location r:=1:                     # radius

 > y3:=gy2(x1,y1,z1,x2,y2,z2,x3,r); z3:=max(gz2(x1,y1,z1,x2,y2,z2,x3,r));

Viewing the Plot of Spheres

 > c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): plots[display]([c1,c2,c3], scaling=constrained,                style=patch, axes=boxed);

The above solutions fail with a zero divide when y1=y2.  Therefore, a new set of equations is derived for this case.

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x2)^2+(y-y1)^2+(z-z2)^2=(2*r)^2: eqn3:=x=x3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}): assign(solution): gy3:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r): gz3:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r): gy4:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gy3(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end: gz4:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gz3(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end:

Another set is needed for the special case of y1=y2 and z1=z2.

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x2)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn3:=x=x3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}): assign(solution): gy5:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r): gz5:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r): gy6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gy5(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi; end: gz6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;   try   res:=[allvalues(gz5(x1,y1,z1,x2,y2,z2,x3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end:

Geometric Model 3 (2 Spheres and Fixed y to Support New Sphere)

Given two spheres located in 3D and fixed y origin location for a third sphere, find the origin of the third sphere that rests on top of the two spheres.  Return FAIL for the z coordinate if it is not feasible.

Let:  r = radius of the spheres

(x1,y1,z1), (x2,y2,z2) = origins of the two given spheres
y3 = known y location of the resting sphere

Then, the origin of the resting sphere is the point where the two given spheres would intersect if their diameters were doubled.  The algebra is solved by Maple as follows:

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2: eqn3:=y=y3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions.  They are converted to Maple functions as follows:

 > assign(solution): hx:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r): hz:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r):

The functions may return multiple values.  Also check for an infeasible solution:

 > hx2:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hx(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end: hz2:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hz(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;    end:

Testing the results,

 > x1:=1.: y1:=1.3: z1:=0:   # sphere 1 x2:=2.5: y2:=1.5: z2:=0:  # sphere 2 y3:=1:                    # sphere 3 y location r:=1:                     # radius

 > x3:=hx2(x1,y1,z1,x2,y2,z2,y3,r); z3:=max(hz2(x1,y1,z1,x2,y2,z2,y3,r));

Viewing the Plot of Spheres

 > c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): plots[display]([c1,c2,c3], scaling=constrained,                style=patch, axes=boxed);

The above solutions fail with a zero divide when x1=x2.  Therefore, a new set of equations is derived for this case.

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x1)^2+(y-y2)^2+(z-z2)^2=(2*r)^2: eqn3:=y=y3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}): assign(solution): hx3:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r): hz3:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r): hx4:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hx3(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end: hz4:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hz3(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;    end:

Another set is needed for the special case of y1=y2 and z1=z2.

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2: eqn2:=(x-x1)^2+(y-y2)^2+(z-z1)^2=(2*r)^2: eqn3:=y=y3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}): assign(solution): hx5:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r): hz5:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r): hx6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hx5(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;   end: hz6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;   try   res:=[allvalues(hz5(x1,y1,z1,x2,y2,z2,y3,r))]   catch: res:=FAIL;   end try;   if res=FAIL then RETURN(FAIL) fi;   if nops(res)=2 and Im(res[1])<>0 then      RETURN(FAIL);   else RETURN(op(res)) fi;    end:

Geometric Model 4 (1 Sphere and Fixed (x,y) to Support New Sphere)

Given fixed (x,y) origin locations and the origin of a sphere, find the origin of the new sphere resting on it.  Return FAIL for the z coordinate if it is not feasible.

Let:  r = radius of the spheres

(x1,y1,z1) = origin of the given spheres
(x2,y2) = known (x,y) location of the resting sphere

The z-coordinate of the origin is found by Maple as follows:

 > unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn:=(x2-x1)^2+(y2-y1)^2+(z-z1)^2=(2*r)^2: solution:=solve(eqn,{z});

A Maple function is programmed for the calculation.

 > iz:=proc (x1,y1,z1,x2,y2,r) local res;   res:=-x2^2+2*x1*x2-x1^2-y2^2+2*y1*y2-y1^2+4*r^2;   if res>0 then RETURN(z1+sqrt(res))      else RETURN(FAIL) fi; end:

Testing the results,

 > x1:=0: y1:=0: z1:=0:   # sphere origin x2:=0.5: y2:=0.5:      # fixed (x,y) of second sphere r:=1:                  # radius of spheres z2:=iz(x1,y1,z1,x2,y2,r);

Viewing the Plot of Spheres

 > c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): plots[display]([c1,c2], scaling=constrained,                style=patch, axes=boxed);

Geometric Model 5 (Distance Between Two Spheres)

Given two spheres in 3D, find the distance between the origins.

Let  (x1,y1,z1), (x2,y2,z2) = origins of the two given spheres .

 > Distance:=(x1,y1,z1,x2,y2,z2)->sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2);

Miscellaneous Routines

Random number generator

 > randomize(): randomX:=stats[random,uniform[0,1]]('generator'):

Box Module

Module Box()  Programming

The Box module performs the calculations of randomly adding balls to the box, keeps track on the balls within the box and provides a means of plotting the results.

Definitions:

Box:-width = width of the box

Box:-depth = depth of the box

Box:-height = height of the box

Box:-radius = radius of the spheres to be placed into the box

Box:-SphereSet = set of data for the spheres within the box

Box:-Spheres = number of spheres within the box

Box:-fxbox, Box:-fybox, Box:-fzbox = routines that link with the calculations of section 'Geometric Model 1' for  a sphere resting on three spheres
Box:-gybox, Box:-gzbox = routines that link with the calculations of section 'Geometric Model 2' for a sphere resting on two spheres and a wall (fixed x)

Box:-hxbox, Box:-hzbox = routines that link with the calculations of section 'Geometric Model 3' for a sphere resting on two spheres and a wall (fixed y)

Box:-makeSphere = routine for adding a sphere to the set.  It tests to make sure that the sphere location is valid.

Box:-checkDistances = routine for checking the distances between a sphere and other spheres.  It must be greater than 2*radius.

Box:-checkDistances2 = another distance checking routine

Box:-firstLayer = routine for randomly placing the first layer of spheres within the box.

Box:-moreSpheres = routine for randomly adding spheres to the box.  It randomly picks three spheres from the set of spheres already within the box and tries to place a sphere on them.

Applying the Stochastic Model

The Box module and supporting functions are used to model the random filling of the box.

First Layer

The spheres of the first layer are placed into the box by random selection of (x,y) coordinates.  It takes many tries to find locations that are not already controlled by existing spheres.  Every 10th try, a random sphere is randomly moved toward the x-axis to make space for more spheres.

 > Box:-firstLayer(5000);

Fill the Box

The box is filled by placing more spheres on top of the existing spheres.  Three spheres are randomly selected from the box and tested to see if sphere can be supported on top of them.  On every 3rd time, tests are performed to see if a sphere could rest upon two of the random spheres and a wall.  On every 10th time, test are performed to see if a sphere could rest on the random sphere and one of the corners.

 > Box:-moreSpheres(40,2000);

`2, S3, S6, S8`

`2, S5, S6, S9`

`9`

`2, S4, S7, S10`

`10`

`2, S10, S3, S11`

`11`

`2, S2, S1, S12`

`12`

`"corner3", S10, S13`

`"corner4", S10, S14`

`14`

`2, S14, S8, S15`

`15`

`2, S9, S12, S16`

`16`

`2, S14, S15, S17`

`17`

`2, S17, S13, S18`

`18`

`2, S12, S11, S19`

`19`

`2, S19, S16, S20`

`20`

`"corner2", S16, S21`

`21`

`2, S9, S15, S22`

`22`

`2, S17, S22, S23`

`23`

`"corner1", S20, S24`

`24`

`2, S19, S11, S25`

`25`

`2, S25, S18, S26`

`26`

`2, S23, S21, S27`

`27`

`2, S26, S20, S28`

`28`

Plot the Spheres Within the Box

Plotting all spheres within the box,

 > Box:-plot();

 >