Dynamical Chaos in a Simple Model of a Knuckleball
DDynamical Chaos in a Simple Model of a Knuckleball
Nicholas J. Nelson, Eric Strauss
Department of PhysicsCalifornia State University, Chico400 W. 1st St.Chico, CA 95929
Abstract
The knuckleball is perhaps the most enigmatic pitch in baseball. Relying on the presenceof raised seams on the surface of the ball to create asymmetric flow, a knuckleball’strajectory has proven very challenging to predict compared to other baseball pitches, suchas fastballs or curveballs. Previous experimental tracking of large numbers of knuckleballshas shown that they can move in essentially any direction relative to what would beexpected from a drag-only trajectory. This has led to speculation that knuckleballsexhibit chaotic motion. Here we develop a relatively simple model of a knuckleball thatincludes quadratic drag and lift from asymmetric flow which is taken from experimentalmeasurements of slowly rotating baseballs. Our models can indeed exhibit dynamicalchaos as long In contrast, models that omit torques on the ball in flight do not showchaotic behavior. Uncertainties in the phase space position of the knuckleball are shownto grow by factors as large as 10 over the flight of the ball from the pitcher to homeplate. We quantify the impact of our model parameters on the chaos realized in ourmodels, specifically showing that maximum Lyapunov exponent is roughly proportionalto the square root of the effective lever arm of the torque, and also roughly proportionalto the initial velocity of the pitch. We demonstrate the existence of bifurcations that canproduce changes in the location of the ball when it reaches the plate of as much as 1.2 mfor specific initial conditions similar to those used by professional knuckleball pitchers.As we introduce additional complexity in the form of more faithful representations of theempirical asymmetry force measurements, we find that a larger fraction of the possibleinitial conditions result in dynamical chaos. The chaotic behavior seen in our simplifiedmodel combined with the additional complex and nonlinear effects likely present in realknuckleballs provide strong evidence that knuckleballs are in fact chaotic, and that thechaotic uncertainty is likely a significant factor in the unpredictability of the pitch. Keywords:
Kunckleball Baseball Dynamical chaos
Email addresses: [email protected] (Nicholas J. Nelson), [email protected] (Eric Strauss)
Preprint submitted to Applied Mathematics and Computation September 14, 2020 a r X i v : . [ phy s i c s . pop - ph ] S e p . “Throwing a butterfly with hiccups” The knuckleball is one of the most specialized pitches in baseball. Most pitcherstry to confuse batters by throwing a variety of pitches such as fastballs, sliders, andcurveballs, each with a unique trajectory achieved by changing the speed and rotationof the ball. Knuckleballers, however, achieve confusion by repeatedly throwing the samepitch. The knuckleball, with its extremely slow rotation in flight, produces a trajectorythat is difficult for pitchers, catchers, and hitters alike to predict. Major League BaseballHall-of-Famer Willie Stargel once stated that “throwing a knuckleball for a strike is likethrowing a butterfly with hiccups across the street into your neighbor’s mailbox” [1].Knuckleballs are thrown in such a way as to minimize the rotation of the ball once ithas left the pitcher’s hand. In the words of long-time knuckleball pitcher Tim Wakefield,“If it spins at all it basically doesn’t work.” [2] Baseballs have a pattern of raised seams,which break the symmetry of the flow of air around the ball and result in aerodynamicforces the depend upon the orientation of the ball relative to its direction of motion.As the ball slowly spins these forces change in both magnitude and direction, producingrelatively large changes in the ball’s vertical and horizontal velocities as it travels. Thismakes knuckleballs very difficult to hit. Similar behaviors are seen in other sports, suchas volleyball and soccer, though in those cases the asymmetric flow is likely related tohydrodynamic instabilities in the flow around the ball.Knuckleballs have a reputation in baseball as being sensitive to how they are throwndown to levels beyond the pitchers ability to control. “You can throw two knuckleballswith the identical release, the identical motion, in the identical place, and one might goone way and the second might go another way” is how professional knuckleball pitcherR.A. Dickey describes them [3]. Due to the complex motions observed in individualpitches and perhaps to the fact that everyone involved with the pitch is unable to pre-dict its motion, some have speculated that kunckleballs may exhibit dynamical chaos.Dynamical chaos is seen in many physical systems [e.g., 4, 5, 6, 7, 8] and is defined bya strong sensitivity to initial conditions. More precisely, dynamical chaos is seen when asystem shows exponential divergence in time for two sets of nearly identical initial con-ditions. Other hallmarks include bifurcations where the changing of a parameter by aninfinitesimal amount results in qualitatively different behavior, and intermittency, wherechaotic regions in parameter space are interspersed with regions that display non-chaoticbehaviors. In this paper we will show that indeed knuckleballs do show clear evidence ofchaotic behaviors, even for our simplified models, as long as those models admit variationin the rotation rate of the ball.In this paper we will develop a set of simple models for the trajectory of a knuckleballas it travels from the pitcher’s hand to home plate. We will then explore solutions ofthese models and show that they yield dynamical chaos as long as torques on the ballare permitted. Chaotic behavior is realized when only one axis of rotation is considered,as well as the more realistic case when two axes of rotation are allowed. The time-scalefor the growth of uncertainty will be shown to depend on both the effective lever armof the torque, and the speed of the pitch. We will further show that simply by alteringthe initial angle and angular velocity of the ball by very small amounts, a pitch withidentical initial position and velocity can arrive at the plate at a wide range of locations.2 . Developing a Knuckleball Model
In this work we follow the common convention used by Major League Baseball’s Stat-Cast system, which sets a coordinate system where home plate in located at the origin[9]. The +ˆ y direction points from home plate to the pitcher’s mound, the +ˆ z direction isupward, and the +ˆ x direction points to the pitcher’s left (batter’s right). The strike zoneis a rectangular portion of the x - z plane that is bounded on the horizontal edges by thesides of home plate located at x = ± .
216 m. The vertical limits of the strike zone aredefined based on the positions of the shoulders, waist, and knees of the batter as theyprepare to swing. Here we choose to adopt a strike zone that extends from z = 0 .
53 mto z = 1 .
19 m, based on long-term averages of called strikes in Major League Baseballgames [10].The trajectory of an object moving through a fluid (such as air) is a classic problemencountered by most physics students in their first physics class. To determine thetrajectory taken, one must identify and quantify all of the relevant forces on the objectand then integrate in time. The aerodynamic forces on a baseball have been studiedextensively both in laboratory settings and using tracking systems deployed in MajorLeague Baseball games [11, 12, 13, 14, 15]. In addition to gravity, most baseballs in flightexperience significant drag forces as well as Magnus or lift forces due to the rotation ofthe ball. Most thrown or batted balls show rotation rates on the order of 100 rad/s[15]. Knuckleballs, however, rotate less than about 20 rad/s. In this regime the Magnusforce is extremely weak and can be ignored, however an additional force arises due tothe changes in flow around the ball caused by the relative locations of the raised seams.We will call this the asymmetry force. Thus the acceleration of a knuckleball can simplybe written as (cid:126)a = (cid:126)F G + (cid:126)F D + (cid:126)F A m , (1)where m is the mass of the baseball. The gravitational force is the simplest of the threewith (cid:126)F G = − mg ˆ z and g the standard 9.8 m/s . The drag force is only slightly morecomplex with a functional form given by (cid:126)F D = − ρπR C D v(cid:126)v , (2)where ρ is the density of the air, R is the radius of the baseball, v is the magnitude of thevelocity, (cid:126)v is the velocity vector, and C D is a dimensionless coefficient that depends onthe details of the flow around the ball. This drag coefficient has been extensively studiedfor baseballs in both laboratory settings and using data from professional baseball games[16, 14, 15]. While there is some weak dependence on velocity for the range of speedscommon in pitches, we have chosen to use a fixed value of C D = 0 .
346 [14]. Our resultsare qualitatively unchanged using values of C D between 0.32 and 0.37.The force due the asymmetric flow of air around the ball caused by the location of theseams is perhaps the most challenging. Functionally it takes the form (cid:126)F A = 12 ρπR C A v ˆ θ × ˆ v , (3)where C A is a non-dimensional coefficient and ˆ θ is a unit vector fixed on a reference pointon the ball. The directional dependence can be simplified using the approximation that3 igure 1: Schematic diagram of the forces acting on a knuckleball moving directly out of the page. Thedrag force (not shown) would be directed purely into the page. Also noted are the angles ψ and φ whichindicate rotation about the ˆ x and ˆ z axes, respectively. the velocity vector is largely in the ˆ y direction, so ˆ θ × ˆ v must have components in only theˆ x and ˆ z directions. The functional form of Eqn. 3 has been verified in a variety of settings[17, 18, 19, 13, 20, 21, 22, 23], however measurement of C A has proven quite difficult. C A depends on the orientation of the ball (cid:126)θ , the angular velocity vector (cid:126)ω , and velocityvector (cid:126)v . The best available data comes from Morrissey [18], who conducted detailedmeasurements using wind tunnel experiments for non-rotating and slowly rotating ballsat a fixed velocity over the full rotation of the ball for two orientations – the two-seamorientation and the four-seam orientation. The knuckleball is commonly thrown so thatthe the two-seam orientation rotates around the x-axis and the four-seam orientationrotates about the z-axis [1].It is likely that the vertical and horizontal components of the asymmetry force bothdepend on the full angle vector, however as a simplification that will allow us to use theexperimental data, we have chosen to consider models where the non-constant terms inthe asymmetry force take the form C A ˆ θ × ˆ v = G ( φ )ˆ x + H ( ψ )ˆ z , (4)where φ represents the angle about the ball’s ˆ z axis, and ψ similarly about the ball’s ˆ x axis. This is a clear simplification as both G and H likely depend on the full angle vector (cid:126)θ , as well as on (cid:126)ω and (cid:126)v . Figure 1 shows a schematic diagram of the coordinate systemand forces discussed above.If we combine Eqns. 1, 2, and 3 and separate our equations into components, we areleft with dudt = − K D u | (cid:126)v | + K A | (cid:126)v | G ( φ ) (5) dvdt = − K D v | (cid:126)v | (6)4 igure 2: The asymmetry force coefficient C A as a function of ball rotation angle in the two-seam (left)and four-seam (right) orientation. Shown are experimental data from Borg & Morrissey [21] (red line)and the filtered data (blue line) which is used in our Type C simulations. The filtered values show goodoverall agreement with the empirical data, while avoiding high wavenumber variations that can causenumerical artifacts in the differential equation solver. dwdt = − K D w | (cid:126)v | + K A | (cid:126)v | H ( ψ ) − g (7)Here we have combined constants from previous equations so that only three remain K D , K A , and g , and defined the velocity vector (cid:126)v = [ u, v, w ]. One obvious simplificationof these equations can be achieved by setting φ = 0. In this case, the dynamics of theproblem will be confined to a single plane. This two-dimensional model will be confinedto the y - z plane if u = 0 initially. In this way we can explore the effect of the asymmetryterm on only one direction of motion at a time. Our treatment of (cid:126)F A has resulted inthe introduction of two new quantities, namely the angles φ and ψ . For the case of acompletely non-rotating ball, these are simply constants, however knuckleballs generallydo show slow rates of rotation and experimental data has shown that the asymmetry forceyields significant torque on the ball in flight [21]. The drag force can also, in principle,exert a torque as well, however we will not consider it in this paper as the torque dueto the drag force has not been experimentally considered independently from the torquedue to the asymmetry force.In this paper we will consider three types of models, which we label Type A, Type B,and Type C. All three types solve Equations 5, 6, and 7, but they differ in their treatmentof the rotation of the ball. Our Type A models assume constant angular velocities abouteach axis, so that (cid:126)θ = (cid:126)ωt + (cid:126)θ . Type A models essentially assume that any torque on theball is insufficient to significantly alter the rotation of the ball in flight.The constant rotation models, however, ignore the torque associated with the forcegenerated on the ball by the asymmetrical air flow around it. A force in the +ˆ x directionshould, for example, produce a corresponding torque on the ball in the ˆ z direction unlessit is acting precisely at the center of mass of the ball. If this force has any effective leverarm it would then produce a torque, leading to a change in the rotation rate. To datethere is no available data on the changes in rotation rate experienced by a knuckleball inflight, so we will introduce an effective lever arm (cid:126)(cid:96) A for (cid:126)F A , such that the torque applied5ariable Value (Physical) Value (Non-dimensional) x y z u v -35.97 m/s -1.07 w φ π rad 0 to 1 ψ π rad 0 to 1 ω z − . . ω x − . . Table 1: Initial conditions used for all simulations unless otherwise noted. to the ball is simply given by (cid:126)τ A = (cid:126)(cid:96) A × (cid:126)F A = (cid:96) A F Az ˆ x − (cid:96) A F Ax ˆ z . (8)We chose to parameterize our effective lever arm in terms of the radius of the ball anda parameter α , defined such that − < α <
1, leading to (cid:96) A = αr . Thus a choice of α = ± α = 0 wouldyield no torque. The evolution equations for the rotation of the ball about the x and z axes, respectively, are d ψdt = − K ψ α | (cid:126)v | G ( ψ ) (9)and d φdt = K φ α | (cid:126)v | H ( φ ) . (10)Here again we have collected all numerical constants into K ψ and K φ .For Type A and B models, we choose a relatively simply functional form for ourangular dependence where G ( ψ ) = C G sin 4 ψ and H ( φ ) = C H sin 4 φ . Figure 2 shows theangular dependence of the asymmetry force, effectively showing C A for two orientationsof a slowly rotating baseball. Examination of the data in Figure 2 shows that bothorientations produce variations that are roughly four-fold periodic over a rotation, andindeed a Fourier transform confirms that the four-fold periodic mode has the most powerin both data sets. For Type A and B models we choose C G = C H = 0 . G ( ψ ) and H ( φ ). Figure 2 shows the experimental data, as well as the smoothed, filtered version ofthe same data that we employ in our models. We have interpolated the experimental dataonto a regular grid in angle and then applied a Gaussian smoothing to avoid small-scalefeatures which will likely produce unreliable numerical results. Additionally, we havesubtracted the mean value from both datasets as a mean value would imply a preferreddirection that would follow the coordinate system under a rotation of the axes. Thesedata are then used to generate a cubic interpolating polynomial to provide an empiricalfunction for our differential equation solver.In summary, we have three groups of models to explore. They are: • Type A: Constant rotation 6 igure 3: Calculated physical distance between models with identical initial conditions but varyinglevels of relative error per timestep. Distances are calculated compared to a model with a relative errorof 10 − per timestep. This figure shows that in order to assure a maximum accumulated numericalerror of less than one part in 10 , we must employ a relative error per timestep of no more than 10 − . • Type B: Torque feedback with sin 4 ψ asymmetry force • Type C: Torque feedback with empirical asymmetry forceFor each group of models, we can either consider 2D motion with rotation about a singleaxis, or 3D motion with rotation about two axes.To solve these equations we employ the SciPy differential equation solving library [24].Specifically, we have use an adaptive 4th/5th order Runge-Kutta solver [25] to solveEqns. 5, 6, 7, 9, and 10. This implementation includes the ability to provide a maximumlevel of relative numerical error per timestep. Initial conditions were selected to berealistic for a knuckleball thrown by a professional pitcher, and are given in Table 1. Forthe initial orientation of the ball ( ψ and φ ) and the initial angular velocity components( ω x and ω z ) a range of values are considered.
3. Measuring Chaos in Our Models
Our three models for the trajectory of a knuckleball are prime candidates for theinvestigation of dynamical chaos as they are strongly non-linear in both the drag andasymmetry forces. These are deterministic equations, so if we provide suitable initialconditions the trajectory of the ball should be calculable for all future time given sufficientnumerical accuracy in our method of solution, and indeed numerical experiments showthis to be the case. Systems like our knuckleball models, however, can exhibit dynamicalchaos and thus show extreme sensitivity to their initial conditions. In our case, thatcould mean that changing the initial orientation or angular velocity of the ball by aninfinitesimal amount may yield large, qualitative differences in outcome when the ballreaches the plate. It is also important to note here that unlike many prototypical chaoticsystems like the damped driven pendulum or the double pendulum, a knuckleball does7 igure 4: Trajectory of a Type A knuckleball with initial conditions ψ (0) = 0 .
106 rad, φ (0) = 0 . ω x (0) = − .
506 rad/s, and ω z (0) = − .
001 rad/s, shown as a function of y position (left panel, redline) as well as snapshots of the location on the x - z plane from the batter’s viewpoint (right panel, stars)where the color of the marker indicates the distance from the plate ranging from the pitcher (light blue)to arrival at the plate (pink). Also shown for comparison is the trajectory of the same pitch withoutthe effects of the asymmetry force (blue line on left, circle markers on right), and the location of thestrike zone (dotted line). Thus even the simplest model considered can show behaviors that qualitativelymatch those seen in real knuckleballs. not exhibit periodic motion. This limits the conceptual tools available to measure chaoticbehaviors.In order to search for chaos, we need to consider the phase space of these models. TypeA models have a three-dimensional phase space, while Type B and C models have a seven-dimensional phase space. To examine the behavior of these trajectories through phasespace, it is necessary to scale our physical quantities and hence non-dimensionalize them.Here we choose the initial y position of the ball as our length scale, the approximate traveltime of the ball to home plate of 0.5 s as our timescale, and 2 π as our angular scale. Wewill use these scale factors in all discussions of phase space distance. Table 1 also showsthe non-dimensional values for our initial conditions. We will denote non-dimensionalquantities with tildes.In order to investigate the presence of dynamical chaos in these models, we need tocalculate the distance in phase space between solutions with differing initial conditions.The phase space vector for the most general case is given by (cid:126)χ = (cid:104) ˜ u, ˜ v, ˜ w, ˜ ψ, ˜ ω x , ˜ φ, ˜ ω z (cid:105) .We can define the distance in phase space between two models with different initialconditions as the Euclidean distance (e.g., distance d is defined by d = ∆ x + ∆ y foran x - y plane) in that space, which we denote ∆ χ . To study potential chaotic behavior,we must be able to evolve two sets of trajectories through phase space which begin veryclose to each other. We have chosen initial perturbations on the order of 10 − .When searching for chaos in a system with a high-dimensional phase space and non-periodic motion through that phase space, the best mathematical tool for finding andquantifying chaos is the Lyapunov exponent. Lyapunov exponents measure the rateat which two trajectories through phase space converge or diverge as they are evolvedthrough time. Dynamical systems such as our models for the knuckleball can be describedby a set of Lyapunov exponents, which describe the local convergence or divergence of a8et of trajectories through phase space in each dimensional of that phase space. Chaoticsystems are those in which at least one of the Lyapunov exponents of the system is positivefor a given choice of parameters [6, 7]. In practice, the largest Lyapunov exponent can bemeasured by determining the exponential growth rate of ∆ χ over a given time interval.Lyapunov exponents have units of s − , and thus represent the number of factors of e thattwo slightly different trajectories in phase space will diverge or converge by per second.Perhaps more practically, they indicate the intrinsic growth of uncertainty in the lo-cation of the model in phase space. As a simple example, if we assume that the phasespace location of a ball leaving the pitcher’s hand could be measured to 0.1% accuracy ineach component of the ball’s velocity, orientation, and angular velocity and that choiceof initial conditions yielded a set of Lyapunov exponents which all had values of 10 s − ,then when the ball crossed home plate the uncertainty may have grown to roughly 15%in each component. In practice, the uncertainty on all components will not grow equally,but as we will show it is not unreasonable to obtain many e -fold increases in phase spacedivergence over the course of a chaotic pitch.It is also important to note that in almost all chaotic systems there are physical limitsbeyond which two systems cannot diverge. These limits are often related to conservedquantities such as mass, momentum, or energy. Even in non-conservative systems such asthe knuckleball, the change in these quantities is limited by the magnitude of the forces,torques, and/or work on the system. This can be referred to as nonlinear saturation ofthe exponential growth, as it is reminiscent of the nonlinear saturation of instabilities influid mechanics.To investigate the exponential growth of uncertainty in these models, we will need tocontrol numerical error so that it is always much smaller than the phase space separationbetween two models we wish to compare. In order to access the exponential separationTo investigate the accuracy of our numerical solver, we conducted a convergence test,comparing the time evolution of a single set of initial conditions at various levels ofrelative error per timestep. Figure 3 shows the physical space distance between a referencesolution computed with a relative error per timestep of 10 − . Given this evolution ofnumerical error, we choose to conduct all models with a relative numerical error of 10 − per timestep and use initial phase space distances between pairs of initial conditions of10 − . This assures that our numerical error will remain much smaller than the phasespace distances we desire to measure.
4. Examples of Knuckleball Trajectories
Having developed a general framework for our three types of knuckleball models, letus first explore if the trajectories they produce show qualitative agreement with observedproperties of the knuckleball – namely that they show substantial movement in flight andthat they work best for very low levels of initial angular velocity.We will first consider the behavior of models with constant rotation, which we labelType A. For these models the angular velocity vector is constant, resulting in a periodicforce on the ball as it travels towards home plate. Though the amplitude of the periodicforcing is proportional to the velocity squared, the dominant component of velocity inthe y -direction means that the periodic movement of the ball in the x - or z -directions haslittle effect on the magnitude of the velocity and hence the amplitude of the forcing. Thesemodels do, however, show substantial effects similar to those seen in real knuckleballs.9 igure 5: Two-dimensional Histograms of the location of simulated knuckleballs with the initial condi-tions given in Table 1 at home plate for Type A (upper left), Type B (upper right), and Type C models(bottom row). Also shown is the strike zone (black dashed line) for reference. For Type C we haveshown the same view as for Types A and B (lower left), as well as a zoomed-out view (lower right) thatshows all simulated plate locations. Types B and C used α = 0 .
5. Initial angles were randomly selectedfrom a uniform distribution between − π/ π/ Figure 4 shows the trajectory of a sample knuckleball with constant rotation, comparedwith the trajectory of a an identical pitch with no asymmetry force. This Type A modelreproduces a key feature of knuckleballs, specifically the ability to change directions inmid-flight. In this example, the knuckleball experiences a strong leftward accelerationinitially which not only stops the ball’s rightward motion but reverses it. The accelerationthen changes sign again, leading to a ball which was moving away from the strike zoneto instead just slip in the lefthand edge. Similarly, this pitch shows roughly only gravita-tional acceleration downward to start, but as the ball approaches the plate its downwardacceleration increases. Thus what looks initially like a pitch in the upper center of thestrike zone changes through the course of its flight to look as if it will be high and left ofthe strike zone. Finally, this pitch turns back to the right and accelerates downward and10 igure 6: Height of the ball at the plate as determined by changing the initial angle ψ and the angularvelocity ω x for Type A knuckleball models with constant angular velocity. All pitches have identicalinitial position and velocity. Knuckleballs show up to 60 cm of movement up (red) or down (blue)due to the asymmetry force when slowly rotating. Thus even the most simple model considered showssubstantial variation due to the initial orientation and spin of the ball. ends up catching the lower left corner of the strike zone. It should be re-emphasized thatthis knuckleball was thrown with identical initial position and velocity to the drag-onlypitch. Similar antics are seen by knuckleballs with a variety of initial orientations androtation rates for Type A, B, and C models.Perhaps the most obvious question to examine with the three knuckleball models toto ask what range of locations as the ball crosses the plate can be achieved. In orderto investigate this we used the same initial position and velocity shown in Figure 4 andran 10 simulations for random initial angles between − π/ π/ simulated pitches arriving at homeplate in the x - z plane for Type A (upper left), Type B (upper right), and Type C (lower)models. The histograms have been normalized so that the color represents the fractionof pitches per bin. All three types show the possibility of substantial movement relativeto a drag-only pitch, though the amount of motion possible varies significantly. Type Asimulations show as far as 70 cm of movement in either the x - or z -directions, leadingto box shaped region with a cross in the center. Type B models show a similar box andcross pattern, but it is much smaller, with only about 25 cm of movement possible. TypeC models exhibit only the cross and show much more potential motion, with pitchesmoving as much as 1.2 m horizontally and 80 cm vertically. This indicates that for TypeC pitches there is significant interaction between the two axes of motion.To further explore the range of possible behaviors for our Type A, B, and C models,we conducted a systematic exploration for a set of 2D models where motion is confined tothe y - z plane and only rotation about the x -axis is considered. We model the trajectories11 igure 7: Left: Height of the ball at the plate for Type B model pitches the same set of initialconditions used in Figure 6 and with α = 1. A drag-only pitch would arrive at the plate with a height of0.96 m (white), while knuckleballs are simulated to vary between 0.79 m (dark blue) and 1.14 m (darkred). Right: Zoom-in of right hand panel highlighting the same region in parameter space shown inright-hand panel of Fig. 11. The sharp transitions between pitches going high (dark red) and low (darkblue) correspond to the arcs of chaotic behavior seen in Fig. 11. In addition to substantial variationin trajectory with initial orientation and spin, these sharp transitions also hint at the possibility ofbifurcation surfaces in phase space. of pitches with − π/ < ψ (0) < π/ − . < ω x (0) < . ± π/ igure 8: Height of the ball at the plate for Type C model pitches in the two-seam configuration withthe same set of initial conditions described in Fig. 7. A drag-only pitch would arrive at the plate witha height of 0.96 m (bright red), while knuckleballs are simulated to vary between 0.78 m (dark blue)and 1.06 m (dark red). Here many potential bifurcation surfaces are seen, providing hints of chaoticbehavior. ditions but varying the initial orientation and angular velocity for Type C models. Asin previous cases, a pitch with only drag forces and gravity would arrive at the platewith a height of 0.96 m (light red). This may seem to conflict with the wider range ofpositions seen in Figure 5 for Type C models, however this may simply be a statementof the rarity of the extreme values. Additionally the extreme values may be a productof the additional complexity when 3D motion and two-axis rotation is considered.Type C models tend to produce more trajectories lower than the drag-only case thanthey do higher trajectories at the plate. Interestingly, pitch tracking data shows a sim-ilar preference for downward motion [13]. As might be expected given the additionalcomplexity of the Type C models, Figure 8 shows rich, complex structure with manybifurcation surfaces. These sharp changes in pitch behavior illustrate why predicting thepath of the ball might be seemingly impossible, and they provide tantalizing hints ofchaotic behavior.
5. Chaos Depends on Torque
While all three types of knuckleball models yield substantial variation in the locationat which they arrive at the plate, the question of if they are chaotic still remains. Forthis we need to compute Lyapunov exponents and look for models where the distancebetween two nearly identical sets of initial conditions shows exponential growth in time.To test for chaotic behavior, we compute a second trajectory with initial conditionsthat have been perturbed by a very small amount. We apply a random perturbationto the initial angles and angular velocity components such that the initial phase spaceseparation between the two initial conditions was 10 − and then tracked the separationbetween the two trajectories in phase space. This is repeated for 10 trajectories with13 igure 9: Separation in phase space of two trajectories with nearly-identical initial conditions. Thenon-dimensional phase space distance ∆ χ grows exponentially from its initial value of approximately10 − through six orders of magnitude within 0.4 s, for a maximum Lyapunov exponent of 37.8 s − . Theexponential growth of the phase space separation shows clear evidence of dynamical chaos. initial angles and angular velocities randomly selected from uniform distributions. Theinitial angles are selected from a uniform random distribution between − π/ π/ ψ and φ . The angular velocity components are also selected from a uniformrandom distribution between − . . ω x and ω z . This process isconducted for all three model types.For Type A models none of the 10 trajectory pairs computed show exponential sep-aration. All models begin with ∆ χ (0) = 10 − and at home plate none of the phasespace separations grow by more than a factor of four. Thus Type A models do not showdynamical chaos for the range of parameter space covering baseball pitches.Type B models show hints of chaotic dynamics in the bifurcation surfaces and in factour search for dynamical chaos shows that Type B models can produce chaos. Aftercomputing 10 trajectories for randomly chosen initial angles and angular velocities,we find that most trajectory pairs do not show exponential growth in ∆ χ , howeverapproximately 1.6% do show at lest two orders of magnitude increase in phase spaceseparation by the time they arrive at home plate. Thus while most choices of initialconditions do not yield chaotic behavior for Type B models, a significant fraction do.The largest such separation in phase space is plotted in Figure 9. This trajectory hasinitial angles ψ (0) = 0 .
696 rad and φ (0) = − .
474 rad, and initial angular velocity com-ponents ω x (0) = 3 .
543 rad/s and ω z (0) = − .
413 rad/s. Thus the pitch makes roughlyone-quarter of a revolution about both axes on its way to the plate. A second trajectorywith a random perturbation with a magnitude of order 10 − was also computed, and thephase space separation grew more than six orders of magnitude as the ball travels towardsthe plate. This yields a maximum Lyapunov exponent of 37 . − . This indicates thatthe phase space separation between two trajectories will increase by a factor of e roughlyevery 26 ms for this pitch. Due to the large positive value of the maximum Lyapunovexponent, this model with these initial conditions clearly demonstrates dynamical chaos.14 igure 10: Largest Lyapunov exponent measured as a function of α (left panel) and as a function of v (right panel). Largest Lyapunov exponents are calculated from 10 trajectories with initial conditions π/ ≤ ψ ≤ π/ ≤ ω x ≤ . α and the initial y component of velocity v . Left Panel: Variation in the largest Lyapunov exponentachieved for a parameter sweep when v = 36 m/s and 0 ≤ α ≤ λ ∝ α . (red line). Right Panel: Same as left panel, but here alpha is held constant at 0.5, while v is varied between 15.0 m/s and 42.0 m/s (blue stars). The bestfit power law (red line) shows λ ∝ v . . This shows how the free parameter α and the initial speed ofthe pitch v impact the growth rate for chaotic uncertainty in the phase space location of the ball. Now that we have unambiguously found a chaotic set of parameters and initial con-ditions, we can explore how the behavior is changed by changing the effective lever armof the feedback torque and the initial velocity of the pitch. One might reasonably ex-pect that chaotic behavior will become more common for faster pitches and similarly forgreater values of α .Our lever-arm parameter is of particular interest as it is essentially a free-parameter inour Type B and C models. It is difficult to estimate what values of α might be expected.Physical scenarios can be envisioned where α might range from -1 (all of the transverseforce is applied at the leading edge of the ball) to 0 (the transverse force is applieduniformly to all points on one hemisphere of the ball) to 1 (all of the transverse force isapplied at the trailing edge of the ball). Given the boundary layer separation which drivesthe motion [18] we might reasonably expect that α >
0, and if the transverse force isuniformly applied to the back quarter of the ball then α = 0 .
5. Borg & Morrissey reportthat in wind tunnels data they observe torques of approximately 4 × − N m, whichwould correspond to α ≈ . α that we consider here, Borg &Morrissey do conclude that “it would be difficult to throw a pitch that did not begin torotate under the influence of aerodynamic shear” [21].In examining Equations 9 and 10, the size of the torque is dependent on two parameters, α and v . As we have stated, α is largely unconstrained but is likely of order 1 /
2, while v can be, and, indeed, is intentionally varied by the pitcher. Both of these parametersshould increase the torque on the ball, leading to a larger nonlinear interaction betweenthe ball’s orientation and its spin. Effectively, increasing either α or v is analogousto amplifying the nonlinearity in the problem that leads to the chaos observed in Type15 igure 11: Left: Maximum phase space separation for Type B model pitches with φ = 0, ω z = 0, − π/ ≤ ψ ≤ π/ − . ≤ ω x ≤ . − .Values range from black indicating non-chaotic behavior to red indicating marginal chaos to whiteindicating clearly chaotic dynamics. Chaotic behavior is achieved in a roughly circular band in parameterspace. Right: Zoom-in of right-hand panel highlighting the region π/ ≤ ψ ≤ π/ . ≤ ω x ≤ . B models. To quantify the effect of both α and v on the resulting chaotic dynamics,we selected a range of initial conditions with π/ ≤ ψ ≤ π/ ≤ ω x ≤ . φ and ω z are set to zero. We then set values for α and v and run 10 combinations of initial conditions. For each phase space sweep, we measure the largestLyapunov exponent achieved in that range of initial conditions for specified values of α and v .Figure 10 shows how the largest Lyapunov exponent achieved in the specified region ofinitial condition space varies with α and v . These numerical experiments demonstratethat the largest Lyapunov exponent is proportional to α . , tantalizingly close to √ α ,and to v . . This also demonstrates that for all non-zero values of the effective leverarm and all initial pitch speeds studied (and those used by knuckleball pitchers [13]) thebehavior of our Type B model can be chaotic.Beyond simply identifying chaos, we can also quantify the timescale of the exponentialdivergence in phase space. For pitches to exhibit significant chaotic effects, the Lyapunovtime scale τ L = 1 /λ must be much shorter than what we term the flight time scale τ F ≈ y /v . When examining how the ratio τ F /τ L varies with α , we see that it variesfrom about 1 when α = 0 .
01 to as much as 17.5 when α = 1. If we set an arbitrary linerequiring at least 5 e-foldings (a factor of ≈ ) to be possible during the ball’s flight,then values of α as small as 0.15 attain this level of chaotic divergence. This then arguesthat only particularly tuned distributions of the asymmetry force can avoid substantiallevels of chaotic behavior.If we now turn to the dependence of the largest Lyapunov exponent with v , we findan intriguing result. The nearly linear dependence of the largest Lyapunov exponent on v leads to a roughly constant ratio of about 12 Lyapunov times per flight time, meaningthat the initial phase space difference has time to grow by a factor of e ( ≈ ) on16ts journey to home plate independent of how fast the ball is thrown. This providesan interesting possible explanation for why knuckleballs seems to work quite well eventhough they are generally thrown much slower than other baseball pitches designed totake advantage of aerodynamic effects. As we previously stated, when selecting initial angles and angular velocities at randomwe found that only 1.6% produced exponentially growing separations in phase space. Toexplore range of possible behaviors as a function of the initial conditions, we now turn toa parameter sweep of possible initial angles and angular velocities while holding v and α constant at 36 m/s and 0.5, respectively. In order to reduce the dimensionality of thespace to explore, we chose to use the 2D form of the models. We explore the behaviorof the system for − π/ ≤ ψ ≤ π/ − . ≤ ω x ≤ . π/ ≤ ψ ≤ π/ . ≤ ω x ≤ . ψ → − ψ and ω x → − ω x .Chaotic behavior is confined to three distinct arcs in this cut through phase space.These arcs are clearly separate when studied at higher resolution (see right-hand panel).Along these ridges initial phase space separations on the order of 10 − can grow fiveorders of magnitude or more in less than 0.5 s. All of the largest Lyapunov exponentsshown in Figure 10 originate from initial conditions along one of these ridges. Thelogarithmic scaling of the color table in Figure 11 serves to highlight the narrowness ofthese ridges, with phase space separations dropping by orders of magnitude when movingeven slightly away from the arcs.The right-hand panel of Figure 11 shows a zoom-in on the region of the initial conditionspace that shows some of the largest divergence in possible trajectories. This is also theregion where large variability is seen for initial rotation rates on the order of the “lessthan a half turn of the ball on its way to the plate” [1]. Here changes in initial rotationrate of as little as a tenth of a radian per second – imperceptible differences to evencurrent high-speed photography – can yield a change in the ball’s position at the plate ofas much as 10 cm. Similar behavior is seen in the horizontal motion of the ball when φ and ω z are allowed to vary, meaning that the modeled behavior could yield several bat-widths of motion in both horizontal and vertical directions when imperceptible changesin the rotation rate and/or initial orientation of the ball are admitted.Having found chaotic behavior in our Type B models, we might expect that the addi-tional complexity of the Type C models would also yield chaotic dynamics. Indeed thatdoes turn out to be the case. As in our Type B models, Type C models also show chaoticbehavior for many choices of initial angle and angular velocity. Figure 12 shows themaximum phase space separation ∆ χ for initial orientations and angular velocities, sim-ilar to those represented in Figure 11. Roughly 12% of models show ∆ χ > − , whichwe interpret as clearly chaotic. In comparing these results with those seen for Type Bmodels, we see some clear similarities in the presence and overall shape of chaotic tracks.We also see that Type C models have significant ranges in parameter space where nochaotic behavior is seen. As expected Type C models show even greater chaotic behaviorthan Type B models. 17 igure 12: Similar to Fig. 11, but for Type C models in the 4S configuration using empirical datafor the angular dependence of the asymmetry force. Phase space separations ∆ χ range from a lowerlimits of approximately 10 − (black) due to uncertainties in the empirical data, to as large as ∆ χ (cid:38) .
6. Conclusion
It is perhaps unsurprising, given the professional baseball talent devoted to them,that knuckleballs should exhibit such complex and richly dynamic behavior even for therelatively simply models presented here. This is even more striking when compared to thereliability with which other pitches can be modeled [12]. In this paper we have presentedsimplified models that do not include dynamic coupling between the ball and the fluidflow, but rather treat aerodynamic forces simply as functions of the ball’s velocity andorientation. We have shown that these models can mimic the strong variability seen inthe location of a given pitch at the plate with very small changes in initial orientationand spin. We have also shown that our models that include torques on the ball dueto the asymmetry force display dynamical chaos for certain regions of parameter spacecorresponding to the steepest gradients in plate location. Finally, we have also shown thatthe chaotic dynamics produced in these models is likely to produce significant uncertaintyin the trajectory of the ball for realistic uncertainties in the initial conditions of the pitch.All three of our model types show the ability to produce variations in the locationof the ball at the plate for small changes in the initial orientation and spin of the ball.Those variations are most pronounced for the Type C models where significant feedbackon the angular velocity of the ball can lead to very dramatic changes in the position ofthe ball at the plate in both the vertical and horizontal directions.Our Type B and C models show strong sensitivity to initial conditions, leading toLyapunov times as small as 26 ms. For Type B models, we have shown that the timescalefor the exponential growth in uncertainties is strongly dependent on the effective leverarm of the torque from the asymmetry force, as well as on the velocity of the pitch.Interestingly, the near-linear dependence of the Lyapunov exponent on the initial velocity18 igure 13: Scatter plot showing the locations at the plate for 1000 pitches using the Type C model.All pitches have initial orientations and angular velocities which vary by less than 0.01%. For reference,the average initial conditions are given in Table 1 with φ = − .
050 rad, ψ = 0 .
330 rad, ω x = 0 . ω z = 1 .
948 rad/s. Also shown for reference is the strike zone (black dotted line). The tracksthrough the plate location space show the range of plate locations achievable with practically identicalpitches thanks to the chaotic nature of our Type C models. means that the magnitude of separation that can be achieved is nearly independent ofthe initial velocity. This may explain why the knuckleball is effective at lower speeds.The large Lyapunov exponents measured in our Type B and C models demonstrate thatinitial uncertainties in the orientation and spin of the ball can grow through as much asseven orders of magnitude over the ball’s flight from the pitcher to home plate.For pitchers and batters, this leads to the conclusion that predicting the trajectoryof a knuckleball is not just difficult – it is functionally impossible. Figure 13 shows thepossible locations at which a knuckleball can arrive at home plate using our Type Cmodel. All trajectories start with identical initial position and velocity, and with initialorientations and angular velocities that vary less that 0.01%. This small variation ininitial orientation and angular velocity grows as the pitches travel towards home plate,leading to the distribution shown. This produces a stark reminder that chaotic does notequate to random. The pitch considered cannot end up anywhere and indeed most ofthe available space is not accessible with these initial conditions. Chaotic does, however,mean that very small initial uncertainties can lead to very large uncertainties when theball reaches the plate.That physical knuckleballs exhibit some level of dynamical chaos has been speculatedand is perhaps unsurprising [13]. What may be of greater interest is that large chaoticuncertainties can be realized with simplified models considered here. Adding additionalnonlinearities and degrees of freedom in the models is not likely to reduce the chaoticbehavior. Additionally, the level of uncertainty which can be achieved from chaoticbehaviors shows that not only are knuckleballs chaotic, but that chaos may play an im-portant role in making the pitch work. Exponential divergence on timescales as small as26 ms and resulting movement of as much as 1.25 m means that not only are knuckleballmodels chaotic, but that the inherent chaotic uncertainty is much larger than the uncer-19ainty in measurement achieved with modern tracking systems. In future work we hopeto use the large numbers of knuckleballs tracked by Major League Baseball’s StatCastsystem to confirm that real knuckleballs show chaotic dynamics similar to those seen inour models.
Acknowledgments
We wish to acknowledge David Kagan for introducing us to the field and pointing usin the right direction, and John Borg for sharing his data. ES was funded by a SummerResearch Award from the College of Natural Sciences at CSU, Chico.
ReferencesReferences [1] D. Clark, The Knucklebook: Everything You Need to Know about Baseball’s Strangest Pitch, 1stEdition, Ivan R. Dee, Chicago, USA, 2012.[2] T. Wakefield, T. Massarotti, Knuckler: My Life with Baseball’s Most Confounding Pitch, HoughtonMifflin Harcourt, Boston, 2011.[3] R. A. Dickey, W. Coffey, Wherever I Wind Up, Blue Rider Press, New York, 2012.[4] E. N. Lorenz, Deterministic Nonperiodic Flow, Journal of the Atmospheric Sciences 20 (2) (1963)130–141. doi:10.1175/1520-0469(1963)020<0130:dnf>2.0.co;2 .URL http://journals.ametsoc.org/jas/article-pdf/20/2/130/3406061/1520-0469 [5] M. Perc, Visualizing the attraction of strange attractors, European Journal of Physics 26 (4) (2005)579–587. doi:10.1088/0143-0807/26/4/003 .URL https://iopscience.iop.org/article/10.1088/0143-0807/26/4/003https://iopscience.iop.org/article/10.1088/0143-0807/26/4/003/meta [6] S. H. Strogatz, Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry,and Engineering, 2nd Edition, CRC Press, Boca Raton, Florida, 2015.[7] R. J. Brown, A Modern Introduction to Dynamical Systems, 1st Edition, Oxford University Press,Oxford, UK, 2018.[8] E. G. Nepomuceno, M. Perc, Computational chaos in complex networks, Journal of Complex Net-works 8 (1) (feb 2020). doi:10.1093/comnet/cnz015 .URL https://academic.oup.com/comnet/article/8/1/cnz015/5481207 [9] M. L. Baseball, Baseball Savant (2019).[10] J. Roegele, The Strike Zone During the PITCHf/x Era, The Hardball Times Baseball Annual 2014(2014) 1–12.[11] R. G. Watts, R. Ferrer, The lateral force on a spinning sphere: Aerodynamics of a curveball,American Journal of Physics 55 (1) (1987) 40–44. doi:10.1119/1.14969 .[12] A. M. Nathan, The effect of spin on the flight of a baseball, American Journal of Physics 76 (2)(2008) 119–124. doi:10.1007/978-0-387-46050-5_5 .[13] A. M. Nathan, Analysis of knuckleball trajectories, Procedia Engineering 34 (2012) 116–121. doi:10.1016/j.proeng.2012.04.021 .[14] D. Kagan, A. M. Nathan, Simplified Models for the Drag Coefficient of a Pitched Baseball, ThePhysics Teacher 52 (5) (2014) 278–280. doi:10.1119/1.4872406 .[15] G. J. Escalera Santos, M. A. Aguirre-L´opez, O. D´ıaz-Hern´andez, F. Hueyotl-Zahuantitla, J. Morales-Castillo, F.-J. Almaguer, On the Aerodynamic Forces on a Baseball, With Applications, Frontiersin Applied Mathematics and Statistics 4 (January) (2019) 1–9. doi:10.3389/fams.2018.00066 .[16] R. K. Adair, The Physics of Baseball, 3rd Edition, Harper Perennial, New York, 2002.[17] R. Watts, E. Sawyer, Aerodynamics of a knuckleball, American Journal of Physics 43 (11) (1975)960–964.[18] M. P. Morrissey, The aerodynamics of the knuckleball pitch: An experimental investigation intothe effects that the seam and slow rotation have on a baseball, Ph.D. thesis, Marquette University(2009).URL http://epublications.marquette.edu/theses{_}open/8
19] T. Asai, K. Kamemoto, Flow structure of knuckling effect in footballs, Journal of Fluids andStructures 27 (5-6) (2011) 727–733. doi:10.1016/j.jfluidstructs.2011.03.016 .URL http://dx.doi.org/10.1016/j.jfluidstructs.2011.03.016 [20] H. Higuchi, T. Kiura, Aerodynamics of knuckle ball: Flow-structure interaction problem on apitched baseball without spin, Journal of Fluids and Structures 32 (2012) 65–77. doi:10.1016/j.jfluidstructs.2012.01.004 .URL http://dx.doi.org/10.1016/j.jfluidstructs.2012.01.004 [21] J. P. Borg, M. P. Morrissey, Aerodynamics of the knuckleball pitch: Experimental measurementson slowly rotating baseballs, American Journal of Physics 82 (10) (2014) 921–927. doi:10.1119/1.4885341 .[22] B. D. Texier, C. Cohen, D. Qu´er´e, C. Clanet, Physics of knuckleballs, New Journal of Physics 18 (7)(2016). doi:10.1088/1367-2630/18/7/073027 .[23] M. A. Aguirre-L´opez, O. D´ıaz-Hern´andez, F. J. Almaguer, J. Morales-Castillo, G. J. EscaleraSantos, A phenomenological model for the aerodynamics of the knuckleball, Applied Mathematicsand Computation 311 (2017) 58–65. doi:10.1016/j.amc.2017.05.001 .[24] T. S. Community, SciPy Reference Guide v1.3.0 (2019).[25] J. R. Dormand, P. J. Prince, A family of embedded Runge-Kutta formulae, Journal of Computa-tional and Applied Mathematics 6 (1) (1980) 19–26. doi:10.1016/0771-050X(80)90013-3 ..