A really common problem in an area of mathematics known as packing optimization is sphere packing: this involves determining how many spheres can be fit into a box, another sphere, a cylinder, etc. given the size of the spheres (they are all the same), and the size of the container. Within this area of mathematics, important historical events to note began in 1611, with Kepler’s Conjecture in which Johannes Kepler proposed that close packing of spheres is the densest possible packing of spheres, and yields an optimal packing density of approximately 74.048% of the cube it is packed into, or π/3√2. Dr. Thomas Hales finally proved this Conjecture in 2015, after years of case checking by a computer.  Another notable conjecture is Ulam’s Packing Conjecture that states that there is no three-dimensional shape with a lower optimal packing density than the sphere, as proposed by Stanislav Ulam. It was proved by Yoav Kallus from Princeton University that Ulam’s Packing Conjecture was true for all point-symmetric convex solids, however not for non-point symmetric convex solids.

I thought I would explore random packing in this article, so I’m going to try and come up with an equation to describe random packing with a cube of sidelength 20, with spheres of different radii. The outline of this article is I’m going to use a computer simulation to run a lot of random simulations of putting spheres into a box, and then take the average for each radius, and then do power and exponential regression on those results to determine the best way of modeling the results. I decided to use Mathematica for this.

Here is a picture to show what the code would generate, when the radius is 0.8 and the cube’s side-length is 20 units:

I tested this code for radii of 0.3-1.0, and here are the numbers of spheres that could fit into cubes of side length 20, with varying radii:

To make it easier to understand the data, here is a scatter plot:

Here, it can be seen that the graph is not linear, but is either a negative exponential function, or a negative power function. To find out, I plotted the points (x, ln y), and (ln x, ln y). If the former graph was linear, an exponential model best fits the points, and modelling can be used to find the pattern; if the latter graph is linear, a power model best fits the points. Here is the plot for (x, ln y):

Here is the plot for (ln x, ln y):

The second plot of points fits a linear model much better than the first. Therefore, a power model best fits the function Radius by Average Number of Spheres Packed. I used both manual power regression, and power regression with the calculator to find the function, where y is the average number of spheres packed, and x is the radius of the spheres being packed.

Power regression:

$y = ax^b$

$a = 679.28711, b = -3.0389171$
$y = \lfloor679.28711{x}^{-3.0389171}\rfloor$

In the modeling, I added the floor function to ensure that all values the function returns are integer values, as only integer spheres can be packed.

And we’re done! However, it’s worth noting that this is just a model based on a small sample set of data, and because it’s random sphere packing, it’s highly unpredictable. Some uses of this model include allowing various manufacturers involved with spherical objects, such as golf ball manufacturers, to create the optimal packaging to best fit a certain amount of golf balls  into a box (when the balls are randomly dispensed).

Here is the link to the code for anyone who is curious: https://mathematica.stackexchange.com/questions/11940/randomly-packing-spheres-of-fixed-radius-within-a-cube