Pinball machines are an arcade staple. This classic game invokes interesting questions regarding the optimization of a pinball machine, including ideal launch angle or peg layout to achieve maximum airtime. In this project, I created a pinball simulation to model a pinball bouncing around a pegboard based on a range of customizable initial conditions. I collected data regarding the relationship between launch angle and the time the ball took to reach the bottom of the board, as well as the corresponding region of impact. I then implemented a simple machine-learning technique to find the optimal pegboard to yield a particular airtime versus launch angle graph result.
Pinball Simulation
Pinball Simulation
The first step in exploring the physics of a pinball machine was to set up a running simulation. I created a board with a range of customizable initial conditions. As the pinball bounces off the pegs and the machine's boundaries, the interaction between the ball, pegs, and walls models real world mechanics.
Initial conditions
Initial conditions
Setting board conditions required initializing various global variables. The basic board setup requires board boundaries in the x and y dimensions and the ball's initial position and velocity. The components of the initial velocity depend upon the launch angle, and its magnitude depends on vScale.
These variables give physical boundary setup, initial position, initial velocity, and simulation duration.
In[]:=
xMax=15;yMax=15;launchAng=5Pi/12;initialPosition={xMax/2,yMax};vScale=3;initialVelocity=vScale*{Cos[launchAng],-Sin[launchAng]};
I then defined other attributes of the peg setup, including dimensions, spacing, and buffering between the pegs and the board borders. collideDist is used to detect collisions between the pinball and a peg, which occurs when the distance between the pinball and peg centers is the sum of their radii.
In[]:=
ballRad=0.05;pegRad=0.5;collideDist=pegRad+ballRad;pegLayout={5,7};buffer=2;pegSpacingX=(xMax-2*buffer)/(pegLayout[[1]]-1);pegSpacingY=(yMax-2*buffer)/(pegLayout[[2]]-1);
Both the pegs and the bounds of the board (walls) have some coefficient of restitution, which measures the elasticity of a collision. A perfectly elastic collision, where all kinetic energy is conserved, would have coefficient of restitution 1; realistically, energy is dissipated in these collisions. The total simulation runtime is broken into smaller time increments, and used to determine the number of iterations that the simulation will run for.
Setting technical aspects including coefficients of restitution, total simulation runtime, time increment duration, and number of iterations.
In[]:=
eWall=0.97;ePeg=0.97;incline=Pi/2;g=9.81*Sin[incline];totalTime=20;dt=0.001;nSteps=Round[totalTime/dt];
Finally, pegCoordSet stores a list of peg center coordinates. The pegs are offset in alternating rows based on the dimensions from pegLayout and spaced based on pegSpacingX and pegSpacingY, maintaining the symmetry of the board.
Generates a list of peg coordinates to fit the pegboard, based on the dimensions specified in pegLayout:
In[]:=
pegCoordSet=Join[Flatten[Table[buffer+{(j-1)*pegSpacingX,(i-1)*pegSpacingY},{i,1,pegLayout[[2]],2},{j,pegLayout[[1]]}],1],Flatten[Table[buffer+{(j-1+1/2)pegSpacingX,(i-1)*pegSpacingY},{i,2,pegLayout[[2]],2},{j,pegLayout[[1]]-1}],1]];
pegs stores the Disk representation of pegs centered at the coordinates determined in pegCoordSet:
In[]:=
pegs=Disk[#,pegRad]&/@pegCoordSet;
Simulation functions
Simulation functions
The function step takes an incremental approach to updating the state of the pegboard through handling changes in the ball's position and velocity due to the acceleration of gravity and collisions with the pegs and walls. Using helper functions, step takes the current position, velocity, and the relevant peg coordinate list as parameters, returning an updated position vector, updated velocity vector, and the peg coordinate list. Note that in this simulation, the peg coordinate list does not change; its utility as a parameter is expanded upon in a later section.
stateAdjust updates the position and velocity vectors using projectile motion kinematics. The relevant peg's coordinates are selected by identifying the peg whose coordinates yield the minimum distance from the ball's position, indicated which peg is closest to the ball. d is taken to be the distance between the ball's center and the peg's center.
The helper function StateAdjust makes basic updates, used within step:
In[]:=
stateAdjust[vel_,pos_,d_,pegCoord_]:=Module[{pegCenter},newVel=vel+{0,-g}dt;newPos=pos+newVeldt;pegCenter=Flatten[Select[pegCoord,#==Flatten[Nearest[pegCoord,newPos][[1]],1]&],1];d=Norm[newPos-pegCenter];{newVel,newPos,pegCenter,d}];
Next the checkDist helper function tests for collisions. If d is less than collideDist, then the ball and peg necessarily overlap. That is, a collision has occurred between the ball and the nearest peg. To calculate the reflected velocity, n is taken to be the normalized vector from the ball's center to the peg, pointing in the direction perpendicular to the tangent line at the ball and peg's point of intersection.
The dot product newVel.n gives the component of the ball's velocity normal to the peg's surface. Then, the ball's velocity is updated according to standard vector reflection: Reversing this normal component of velocity, and scaling by the peg's coefficient of restitution.
Updating newPos corrects error by resetting the ball position to the collision distance away from the peg's center. Scaling collideDist by n ensures that the offset is given in the proper direction, so that the ball is properly tangent to the peg. This is necessary since there was some small overlap between the ball and peg regions to satisfy the condition d<collideDist in the first place.
The helper function checkDist calculates the ball's position and velocity using standard vector reflection after a collision with a peg:
In[]:=
checkDist[d_,n_,newPos_,center_,vn_,newVel_]:=If[d<collideDist,n=(newPos-center)/d;vn=newVel.n;{newVel-(1+ePeg)vnn,center+ncollideDist},{newVel,newPos}];
checkBounds checks the edge cases for the ball in relation to the board. If the ball has exceeded the maximum or minimum of either dimension and thus exited the board, its position is corrected so the ball is tangent to the edge of the board and its velocity in the offending dimension is reversed and scaled by the wall's coefficient of restitution.
This helper function checks the ball's position in relation to board edges, making position and velocity adjustments.
In[]:=
checkBounds[{pos_,vel_}]:=Which[pos[[1]]>xMax-ballRad,{{xMax-ballRad,pos[[2]]},{-eWall*vel[[1]],vel[[2]]}},pos[[1]]<ballRad,{{ballRad,pos[[2]]},{-eWall*vel[[1]],vel[[2]]}},pos[[2]]<ballRad,{{pos[[1]],ballRad},{vel[[1]],-eWall*vel[[2]]}},pos[[2]]>yMax-ballRad,{{pos[[1]],yMax-ballRad},{vel[[1]],-eWall*vel[[2]]}},True,{pos,vel}];
Taken together as the sum of its parts, step makes multiple updates to the ball’s position and velocity as it checks through the cases described above.
The step function in its entirety:
In[]:=
step[{pos_,vel_,pegCoord_}]:=Module[{newPos,newVel,n,vn,d,pegCenter},{newVel,newPos,pegCenter,d}=stateAdjust[vel,pos,d,pegCoord];{newVel,newPos}=checkDist[d,n,newPos,pegCenter,vn,newVel];{newPos,newVel}=checkBounds[{newPos,newVel}];{newPos,newVel,pegCoord}];
Next, I repeatedly applied the function step to its own output, starting with step as a function of initial conditions. This process runs step on its own updated state, handling the scope of the simulation. Therefore, trajectory is a list of lists. Each sublist records the position, velocity, and peg coordinate list information for each step, though the peg coordinates remain fixed throughout.
Creates a list of information from each step of the ball's simulated trajectory through the pegboard space.
In[]:=
trajectory=NestList[step,{initialPosition,initialVelocity,pegCoordSet},nSteps];
Takes the positions only from each entry in the trajectory list.
In[]:=
positions=trajectory[[All,1]];
Finally, the ball’s position as a function of time is overlaid onto a graphical representation of the pegboard.
This Manipulate generates an animation of the ball’s motion.
In[]:=
Animate[Show[Graphics[{PointSize[ballRad],Point[positions[[Round[frame]+1]]]},Axes->True,AxesOrigin->{0,0},PlotRange->{{0,xMax},{0,yMax}}],Graphics[{Lighter[Brown],pegs}],ImageSize->600],{frame,0,nSteps,1},SaveDefinitions->True]
Out[]=
Airtime and End Position Versus Launch Angle
Airtime and End Position Versus Launch Angle
Next I plotted the trajectories, airtime, and end x-position for a ball launched over a range of angles given a certain pegboard. This enables visualization and analysis of the consequences of making small adjustments in both peg setup and pinball launch angle.
Airtime versus Launch Angle
Airtime versus Launch Angle
The pinball simulation runs for a certain number of steps, which is directly related to the simulated duration: 1,000 steps is equivalent to one second. Call the “airtime” the number of steps that pass before the ball first hits the bottom of the board, and consider to the game to be terminated at the first occurrence of this event. Airtime depends on the ball’s initial velocity and the board setup, and may be maximized through interactions with board obstacles: For example, the ball falling and bouncing upward off a peg will yield longer airtime than a pure freefall from the starting position.
angTrajPairs gives a list of {launch angle, trajectory, peg coordinate list} triples, where each trajectory is a function of the corresponding launch angle. This essentially repeats the process described for trajectory.
The launch angle ranges from 0 to Pi radians, measured clockwise from the rightward horizontal. The ball’s launch range is controlled by modifying vScale, which may increase the initial horizontal velocity component. The length of angTrajPairs depends on how many sample launch angles are used, determined by the parameter numPaths.
Overall, angTrajPairs stores many trajectories. A single trajectory is a list of lists, because many steps comprise a single trajectory, and each step contains the ball’s position and velocity information after each time increment.
The launch angle ranges from 0 to Pi radians, measured clockwise from the rightward horizontal. The ball’s launch range is controlled by modifying vScale, which may increase the initial horizontal velocity component. The length of angTrajPairs depends on how many sample launch angles are used, determined by the parameter numPaths.
Overall, angTrajPairs stores many trajectories. A single trajectory is a list of lists, because many steps comprise a single trajectory, and each step contains the ball’s position and velocity information after each time increment.
Provides a list organized into discrete sublists for each launch angle, each containing the launch angle, resulting trajectory, and peg coordinate list.
In[]:=
angTrajPairs[numPaths_Integer,pegCoord_List]:=Table[{i,NestWhileList[step,{initialPosition,{vScale*Cos[i],-Sin[i]},pegCoord},#[[1,2]]>ballRad&]},{i,0,Pi,Pi/numPaths}];
Next, I wrote a few functions and values to store information about a pegboard and its trajectories.
Calls angTrajPairs for a given peg coordinate set:
airtimePairs stores modified launch angle and trajectory information. Each element is of the form {launch angle, corresponding trajectory length}.
Extracts the relationship between launch angle and trajectory length:
results takes the second element from each of the items in angTrajPairs, meaning it accesses each individual trajectory in sequential order and excludes the launch angles and peg coordinate lists.
Extracts the list of trajectories:
showTrajectories displays the ball's "trail" by accessing all position coordinates from the list of trajectories and overlaying this result on top of the pegboard graphic to represent the ball's interactions with the board.
Traces out a ball's various trajectories when launched over a range of angles, overlaying the trajectories onto the pegboard of focus.
Displays the trajectories and resulting airtime graphs:
End X-Position Versus Launch Angle
End X-Position Versus Launch Angle
The second graphing focus for this simulation is final x-position, or termination position, versus launch angle. endXPositions accesses the final x-position only from each trajectory, meaning the x-position from the last step in each trajectory.
Final x-positions from each trajectory:
endPairs re-inserts corresponding launch angle to obtain {launch angle, final x-position} pairs.
This information is sufficient to plot the launch angle and final x-position. Note the rotational symmetry about x=Pi/2, which is due to the symmetry of the pegboard and trajectories.
Final x-position versus launch angle graph:
End X-Position Distribution
End X-Position Distribution
Consider the board to be divided into some number of equally-sized horizontal regions. endRegion determines which of these sections the ball landed in. Formatting termination data as a bar graph better represents the distribution of terminations across termination sites.
The numerical region of termination based on ending x-position:
Set the number of regions:
Frequency maps over the lists of final x-positions to determine the total count of terminations per end region.
A tally of how many terminations occurred in each end region:
However, frequency accounts only for regions with a nonzero number of terminations. To adjust for this, zeroCases finds the regions with zero terminations and assigns 0 to each of those values following frequency’s {region number, termination count} pair format.
Creates filler values for regions with zero terminations:
Joins the elements of frequency an zeroCases, ordering numerically by region number:
I displayed this data in a bar chart, with the region number on the x-axis and the number of terminations within that region on the y-axis.
Bar chart recording termination distribution:
Example Outputs
Example Outputs
I applied the steps derived above to a range of peg arrangements, velocity scales, and number of launch angles:
Optimizing a Pegboard to Fit an Airtime Graph
Optimizing a Pegboard to Fit an Airtime Graph
The methods above take a certain pegboard and generate airtime and final x-position graphs. This raises the question of approaching the problem from the other direction: Starting with a pegboard and its respective airtime graph, then incrementally modifying peg position until the pegboard’s airtime graph best fits some predetermined target end graph.
This sort of reverse-engineering may be achieved through some measure of error calculation between graphs and implementing a machine-learning method to reduce that error. My approach focused on local optimization. I calculated the error of a current state, determined the error of multiple possible subsequent states, and moved to the state of greatest relative error reduction. Here, a state change involves moving a peg in one of the four principal orthogonal planar directions. An optimal state is supposed to be reached once peg positioning has resulted in a local error minimum.
This sort of reverse-engineering may be achieved through some measure of error calculation between graphs and implementing a machine-learning method to reduce that error. My approach focused on local optimization. I calculated the error of a current state, determined the error of multiple possible subsequent states, and moved to the state of greatest relative error reduction. Here, a state change involves moving a peg in one of the four principal orthogonal planar directions. An optimal state is supposed to be reached once peg positioning has resulted in a local error minimum.
vScale, numPaths, and pegCoordSet must be customized to test different scenarios.
Below, I used a sine wave as my target function model, and plotted multiple points along the sine wave.
Scaling the output to roughly match the scale of a typical trajectory:
Setting and scaling the target end function:
refPoints generates a set of target reference points, where the x-values match the launch angles used for the trajectory that the target will be compared against.
Mapping the target function onto a finite number of x-values:
The actual {launch angle, airtime} values that the optimization method aims to fit to the sine wave model are generated using angTrajPairs.
Gives a list of actual launch angle versus airtime coordinates for a set pegboard:
Error calculation
Error calculation
I utilized root mean square error (RMSE) as my error calculation method. RMSE is calculated as the the square root of the average squared difference between corresponding values from two functions. In this case, I applied RMSE to pegAirtimeResults and refPoints as the input and target lists, and assigned the result to baselineError.
Computes root mean squared error, a standard metric used in machine learning.
Calculates the baseline error of the original pegboard:
Error minimization
Error minimization
Next, I aimed to minimize the RMSE value at any given step. The function greedyMove implements a greedy algorithm on a list of pegboard coordinates. A greedy algorithm is a decision-making approach focused on choosing the locally optimal solution at each step.
One iteration of greedyMove changes the location of a specified peg on the pegboard, reducing the error between the pegboard’s airtime function and the target function.
One iteration of greedyMove changes the location of a specified peg on the pegboard, reducing the error between the pegboard’s airtime function and the target function.
Helper functions
Helper functions
findMin is a helper function for greedyMove. This function returns the minimum RMSE value from the elements within options, and index within options at which this value was located. options is a list within greedyMove, which is defined below its helper functions. options gives four peg adjustment choices: Moving the peg position one unit in either the positive or negative x or y direction.
minErrorIndex first assumes that the minimum error is found using the first element of options. The function then iterates through the remaining elements of options to determine whether the other possibilities further minimize error, updating minErrorIndex accordingly.
findMin returns the minErrorIndex pair, which contains the minimized error and the index at which this min error is found within options.
minErrorIndex first assumes that the minimum error is found using the first element of options. The function then iterates through the remaining elements of options to determine whether the other possibilities further minimize error, updating minErrorIndex accordingly.
findMin returns the minErrorIndex pair, which contains the minimized error and the index at which this min error is found within options.
Helper function called by greedyMove to determine which peg movement direction minimizes error:
The function testBest compares the new error to the previous state’s error. If the change identified by findMin gives an error lower than the current state’s error, testBest updates the pegboard to that lower-error state. Otherwise, the pegboard maintains its current state. This effectively judges whether the next best minimized error is actually optimal.
Helper function called by greedyMove[] to control updating the pegboard and error value, if applicable.
Optimal adjustments and recursion
Optimal adjustments and recursion
greedyMove takes as parameters the current peg coordinate list, the current error value, and the index within the peg coordinate list indicating which peg will be adjusted. The function returns the updated corresponding values. greedyMove uses findMin and testBest to determine the optimal direction, if any, of peg position adjustment.
Changes the state of the pegboard to the lowest-error state of the four nearest neighboring states.
Applying greedyMove recursively repeatedly adjust peg position to a lower-error state until the RMSE value has reached a local minimum. This minimum indicates that the RMSE of all immediate neighboring states is either higher or equal to that of the current state, and therefore cannot be improved from a local standpoint. At this point the recursion terminates.
Note that the greedy algorithm will not necessarily reach the globally optimal solution. However, it will obtain a satisfactory approximation at a local error minimum.
Note that the greedy algorithm will not necessarily reach the globally optimal solution. However, it will obtain a satisfactory approximation at a local error minimum.
Recursion using greedyMove, starting with initial peg coordinate list and error.
Extracts the final, optimized peg coordinate set from the list of states:
This Manipulate demonstrates the peg adjustment process for a single peg.
Displays the pegboard updating after each iteration of greedyMove:
Then I displayed the target airtime graph side by side with the initial and optimized airtime graphs, to compare the results. For this single-peg case, the optimized airtime graph resembles the sine wave. The two graphs compare results from before and after applying the greedy algorithm.
First, this general helper function displays points from a function depending on launch angle:
Graphics output of the target airtime graph versus the starting actual airtime graph:
Graphics output of the target airtime graph versus the final actual airtime graph:
Conclusions
Conclusions
In this project, I created algorithms to successfully simulate a pinball board, modeled the functional dependency of airtime and game termination location as a result of various initial conditions, and implemented a method to find an optimal pinball board given an airtime graph. Each subsequent piece of this project built upon the previous pieces, producing informative visualizations and leading to board optimization. This project demonstrated that not only is it possible to modify a pegboard to approximately fit a particular airtime graph, but this optimization works decently well on simple and complex pegboards alike.
Future Directions
Future Directions
Limited time prevented me from implementing a more rigorous machine learning or global optimization method. Therefore, steps to further develop this project include training a neural network or using stochastic gradient descent as part of a fine-tuned peg optimization method, finding the best way to iterate through pegs on a multi-peg board, and ultimately reaching a globally optimal solution. I would also like to examine the application of the target-graph fitting process not only to airtime, but to final termination x-position and distribution.
References
References
◼
Ball Bouncing on Hilly Terrain. (2017). Mathematica Stack Exchange. Mathematica Stack Exchange. https://mathematica.stackexchange.com/questions/141209/ball-bouncing-on-hilly-terrain
Acknowledgements
Acknowledgements
First, I would like to thank my advocate and mentor Piero Sanchez for his guidance and insights throughout this project, and for joining our mentor group on countless hacky sack breaks. Thank you also to the TAs, particularly Hendry Xu for offering approaches to optimization; to Stephen Wolfram for helping me define a goal for this project; and to directors Eryn, Megan, Rory, and Cyrus for facilitating the incredible experience that has been WSRP 2026.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Pinball machine mechanics and peg position optimization
by Isabelle Archer
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3752547
by Isabelle Archer
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3752547