Mass Dependent Computational Analysis of Projectile Motion under Quadratic Air Drag Using the Runge-Kutta Method

Annasi Ayubu Said^1*, Haroub Pollycap Mbewe^2*, Melckzedeck Michael Mgimba³, Hassani Saidi Namanolo⁴, Shadya Muhammed Rashid², Sumaiya Issa Ussi²

¹Department of Natural Sciences, College of Sciences and Technical Education, Mbeya University of Science and Technology, Mbeya, Tanzania.
²Department of Natural Sciences, Faculty of Science, Abdulrahman Al-Sumait University, Zanzibar, Tanzania.
³Department of Geosciences and Mining Technology, College of Engineering and Technology, Mbeya University of Science and Technology, Mbeya, Tanzania.
⁴Department of Mathematics and Computer Science, Faculty of Science, Abdulrahman Al-Sumait University, Zanzibar, Tanzania.

DOI: 10.4236/ojapps.2025.1512260                

Abstract

This study investigates projectile motion under quadratic air drag, focusing on mass-dependent dynamics using the Runge-Kutta (RK4) method implemented in FreeMat. Quadratic drag, predominant in high Reynolds number flows, significantly alters trajectories and challenges analytical modeling. We examine projectile trajectories for varying masses: 0.046 kg (golf ball), 0.145 kg (baseball), and 0.189 kg (softball) across initial velocities ranging from 10 m/s to 50 m/s. Simulations show that increasing mass yields greater resistance to deceleration, resulting in longer flight ranges and higher terminal velocities. Distinct asymmetries are observed in vertical velocity profiles during ascent and descent, validating theoretical expectations. The study quantifies the influence of mass on terminal velocity, range, and peak height, providing insights applicable to ballistic modeling, sports physics, and atmospheric transport phenomena. The RK4 method demonstrates accuracy and stability in solving nonlinear differential equations governing such dynamics. These findings emphasize the significance of quadratic drag modeling in realistic projectile simulations and provide a validated numerical framework for applications in ballistics, sports physics, and atmospheric dynamics.

Keywords

Projectile Motion, Air Resistance, Quadratic Drag, Nonlinear Dynamics, Numerical Simulation, Mass Dependence

1. Introduction

Accurate modeling of projectile trajectories in resistive media is essential in disciplines ranging from sports physics to aerospace engineering and environmental science [1]-[5]. While projectile motion in a vacuum is well understood and follows simple parabolic trajectories, the presence of air resistance, particularly when it is proportional to the square of velocity, introduces significant nonlinearities. These nonlinear effects often render closed-form analytical solutions intractable [6]-[8]. Despite centuries of investigation, the dynamics of projectiles in the presence of quadratic air drag remain an open area of applied physics, especially when considering mass-dependent effects [9]-[11].

Quadratic drag becomes dominant when projectiles move at moderate to high Reynolds numbers, where inertial forces exceed viscous forces. This is typical in real-world conditions involving spherical bodies such as baseballs, golf balls, and softballs [12]-[14]. For such scenarios, linear drag models—which assume drag is proportional to velocity—fail to capture the realistic aerodynamic behavior, especially when velocities exceed a few meters per second [13] [15]. In contrast, quadratic drag models more accurately reflect empirical observations and are therefore preferred in simulations involving sports, ballistics, and fluid-mediated motion.

Several studies have approached this problem using approximate analytical methods, such as perturbation theory, truncated power series expansions, and the homotopy analysis method [9] [16]-[18]. While these methods offer partial solutions, they often suffer from poor accuracy outside small perturbation domains or short time intervals. The resulting expressions are typically too complex to allow generalization or physical interpretation [19]. As a result, numerical methods have become essential tools for capturing the full complexity of projectile dynamics under nonlinear drag forces [20]-[22]. Among these, the fourth-order Runge-Kutta (RK4) method is widely regarded for its balance of accuracy and computational efficiency in solving coupled nonlinear differential equations [23].

Nevertheless, a key limitation in existing literature is the underrepresentation of mass sensitivity in projectile motion under quadratic drag. Most prior works either focus on a single mass or treat mass as an implicit scaling parameter rather than an explicit variable influencing the dynamics. This omission limits the application of such models to real-world systems where object mass varies significantly. Understanding how mass affects projectile trajectories under nonlinear drag is particularly crucial for optimizing performance in sports, predicting the motion of airborne particulates, and developing robust aerospace trajectory simulations [6] [11] [24] [25].

This study addresses these limitations by incorporating mass-dependent effects into the modeling of projectile motion under quadratic drag. Using the fourth-order Runge-Kutta method implemented in FreeMat, we numerically solve the nonlinear differential equations governing projectile dynamics for three distinct spherical masses representing a golf ball (0.046 kg) [26], a baseball (0.145 kg) [27], and a softball (0.189 kg) [28]. We systematically analyze how varying mass influences the projectile’s range, flight time, vertical displacement, terminal velocity, and the shape of its trajectory.

The results highlight distinct asymmetries in vertical velocity profiles between ascent and descent, and demonstrate that heavier objects experience slower deceleration and greater horizontal displacement. Our simulations also validate theoretical expectations regarding drag dominance and terminal behavior. This work presents a validated numerical framework for studying mass-sensitive projectile motion under realistic aerodynamic conditions. The findings have broad implications across multiple fields, including sports physics, environmental modeling, and defense applications.

The other part of this study is organized as follows: Section 2 presents the derivation and discussion of equations governing projectile travel under quadratic drag force, taking into account the projectile’s parameters, as well as the horizontal and vertical components prescribed by Newton’s second law of motion. Section 3 presents the analysis and discussion of the obtained results, while the conclusion and summary of this work are presented in Section 4.

2. Methodology

2.1. Quadratic Drag and Mass Dependence

In real-world projectile motion, especially at high speeds, the resistive force exerted by air is more accurately modeled as being proportional to the square of the velocity. This quadratic drag model is well-supported in classical mechanics and is essential for accurate predictions in atmospheric flight, sports physics, and ballistic trajectories.

As established by Marion and Thornton [13] and further elaborated by Kane and Levinson [29], the drag force $F_d$ acting on a body moving through a fluid at high Reynolds number is given in Equation (1) below.

$F_d=\frac{1}{2}{C}_d\rho A{v}^2$ (1)

where $C_d$ is the drag coefficient, $\rho$ the fluid density, $A$ the cross-sectional area, and $v$ the velocity. This quadratic dependence becomes dominant when inertial forces transcend viscous effects, typically at speeds exceeding a few m/s for small projectiles in air.

Lubarda [30] and Carmago [31] provided complete analytical and numerical treatments of projectile motion under quadratic drag, demonstrating its nonlinearity and sensitivity to both initial conditions and mass. Their work aligns with our simulation results showing that lighter projectiles are significantly more affected by air resistance, due to their lower momentum and higher deceleration rates under equal drag forces.

Specifically, since the acceleration due to drag is inversely proportional to mass, as in Equation (2):

$a_d=\frac{F_d}{m}\alpha \frac{v^2}{m}$ (2)

A smaller mass results in a larger deceleration for the same velocity. This explains the reduced vertical displacement observed in our simulations for the 0.046 kg projectile, as compared to the 0.198 kg projectile, over equal time intervals.

This study presents a numerical investigation of vertical and horizontal components under the influence of gravity and quadratic air resistance. The motion is modeled using classical Newtonian mechanics, with particular attention to nonlinear drag forces that significantly affect objects moving at high velocities through a fluid medium. Due to the complexity of the resulting differential equations, a computational approach is adopted [32]. All simulations and numerical computations were performed using FreeMat, a high-level interactive environment for numerical computation and data visualization [33].

The use of FreeMat allowed for efficient coding, execution, and visualization of the results, including time-series plots of velocity and displacement [34] [35]. The numerical accuracy and stability of the simulation were ensured through careful selection of the integration time step and validation against known theoretical behavior in limiting cases (e.g., motion in vacuum or low-speed regimes) [36] [37]. Given the nonlinearity of the differential equations, closed-form analytical solutions are generally not obtainable. Therefore, a numerical solution is implemented using FreeMat. The simulation employed a time-stepping algorithm based on the classical fourth-order Runge-Kutta (RK4) method to approximate the velocity and displacement of the object as functions of time [21] [38].

Key simulation parameters, such as initial velocity, object mass, drag coefficient, and time step, were defined and systematically varied to evaluate their influence on the motion. The computational model tracks the object’s velocity and position over time [39] [40], enabling the analysis of critical phenomena such as terminal velocity and time to peak height. However, Simulations used a fixed time step of Δt = 0.01 s, chosen based on a convergence analysis showing that halving Δt altered the final range by less than 0.1%, confirming solution insensitivity to further refinement and ensuring both accuracy and computational efficiency.

2.2. Aerodynamic Parameters

The aerodynamic parameters governing drag, namely the drag coefficient ($C_d$ ), air density ($\rho$ ) was taken as $\rho =1.225\text{kg}/{\text{m}}^{\text{3}}$ [41], and the cross-sectional area ($A$ ), directly influence the drag constant $c$ , which appears in the quadratic drag formulation in Equation (3) below.

$c=\frac{1}{2}{C}_d\rho A$ (3)

Table 1 summarizes these parameters for the three projectiles examined. The values of $C_d$ were selected from experimentally established ranges for spherical bodies moving at moderate to high Reynolds numbers ($Re~{10}^4\text{-}{10}^5$ ) [42] [43], ensuring that the assumptions of Reynolds numbers independent drag is appropriate for the velocities analyzed. The diameters correspond to regulation sporting equipment and were used to compute the frontal area via $A=\pi d^2/4$ .

Table 1. Aerodynamic parameters used.

Projectile	Mass (kg)	Diameter (m)	Cd	Area A (m2)	c (kg/m)
Golf ball	0.046	0.0427	0.47	1.43 × 10−3	0.000 41
Baseball	0.145	0.073	0.30	4.19 × 10−3	0.000 77
Softball	0.198	0.097	0.40	7.39 × 10−3	0.001 81

The hierarchy of drag constants and $\frac{c}{m}$ ratios in Table 1 is consistent with classical aerodynamic scaling: drag force increases with frontal area and $C_d$ , whereas resistance to deceleration increases with mass [10] [44] [45]. However, heavier projectiles preserve momentum longer and exhibit less curvature distortion and range suppression under quadratic drag.

Moreover, the aerodynamic parameters in Table 1 create the physical foundation for the observed mass-dependent behavior in subsequent results. By accurately representing realistic projectile geometries and aerodynamic loading conditions, these parameters ensure that the simulated trajectories faithfully capture the nonlinear dynamics of projectiles moving through atmospheric flows.

2.3. Horizontal Component

In many real-world scenarios, especially at higher velocities, objects moving through a fluid (such as air or water) experience a drag force that is proportional to the square of their velocity. This type of resistance is referred to as quadratic drag, and it plays a significant role in the dynamics of projectiles, vehicles, and objects in free fall. In this section, we consider the horizontal motion of an object subject to such a drag force and derive analytical expressions for its velocity and displacement as functions of time.

To begin the analysis, we apply Newton’s second law in the horizontal direction, considering only the drag force, which acts to oppose the motion. The force due to quadratic drag is given by $F=-c{v}^2$ , where $v$ is the velocity, and ccc is the drag coefficient, which depends on the shape, surface area, and properties of the medium. However, for a standard size, it is the quadratic drag force that affects the motion. For this reason, the equation of motion (in terms of $v$ ) becomes as in Equation (3).

$m\dot{v}=mg-c{v}^2$ (4)

Equation (4) can also be written with components from vertical and horizontal as.

Since $g=0$ , then

$m{\dot{v}}_x=-c\sqrt{v_x^2+v_y^2{v}_x}$ (5)

$m{\dot{v}}_y=mg-c{v}^2=mg-c\sqrt{v_x^2+v_y^2{v}_y}$ (6)

As we noted earlier, these Equations (5) and (6) are coupled and are generally not solvable analytically (in terms of equations); they can be solved numerically.

Applying Newton’s second law $F=ma$ , we obtain the differential equation governing the motion:

$\frac{\text{d}v}{\text{d}t}=-\frac{c}{m}{v}^2$ (7)

This is a first-order nonlinear differential equation due to the $v^2$ term. Nonlinear equations are generally more complex to solve than linear ones, but this equation is separable, which simplifies the process. By rearranging terms, we separate the variables $v$ and $t$ :

$\frac{\text{d}v}{v^2}=-\frac{c}{m}\text{d}t$ (8)

We then integrate both sides of the equation. The left-hand side is integrated with respect to velocity from the initial velocity $v_0$ to $v\left(t\right)$ , and the right-hand side with respect to time from 0 to $t$ :

${\displaystyle \underset{v_0}{\overset{v\left(t\right)}{\int }}\frac{\text{d}{v}^{\prime }}{{v^{\prime }}^2}}=-\frac{c}{m}{\displaystyle \underset{0}{\overset{t}{\int }}\text{d}{t}^{\prime }}$ (9)

The integration yields:

$\left(\frac{1}{v\left(t\right)}-\frac{1}{v_0}\right)=\frac{ct}{m}$ (10)

Solving for $v\left(t\right)$ , we obtain

$v\left(t\right)=\frac{v_0}{1+\frac{t}{\tau }}$ (11)

Where we define a characteristic time constant $\tau$ as:

$\tau =\frac{m}{c{v}_0}$ (12)

This characteristic time $\tau$ indicates how quickly the velocity decays due to drag. Physically, it represents the time it takes for the object’s velocity to reduce significantly from its value.

Next, to find the position $x\left(t\right)$ of the object as a function of time, we integrate the velocity:

$x\left(t\right)=x_0-{\displaystyle \underset{0}{\overset{t}{\int }}v\left(t^{\prime }\right)\text{d}{t}^{\prime }}=x_0+{\displaystyle \underset{0}{\overset{t}{\int }}\frac{v_0}{1+\frac{t^{\prime }}{\tau }}\text{d}{t}^{\prime }}$

This integral evaluates to:

$x\left(t\right)=x_0+v_0\tau \ln \left(1+\frac{t}{\tau }\right)$ (13)

Assuming the initial position $x_0=0$ , the final expression for displacement becomes:

$x\left(t\right)=v_0\tau \ln \left(1+\frac{t}{\tau }\right)$ (14)

This function, as described in Equations (14) and (11), represents the horizontal velocity and position of an object experiencing quadratic drag as it decelerates over time, respectively. Unlike motion in a vacuum, where velocity remains constant in the absence of horizontal forces, the quadratic drag causes a rapid initial deceleration, which gradually slows as the object continues moving and its velocity decreases. The logarithmic nature of the position function reflects the diminishing displacement over time as drag dominates the motion.

2.4. Vertical Component

In the vertical direction, the motion of an object subjected to both gravitational force and quadratic air resistance is governed by the following equation of motion for downward (15) and upward motion (16):

$m\frac{\text{d}v}{\text{d}t}=-mg-c{v}^2$ (15)

$m\frac{\text{d}v}{\text{d}t}=-mg+c{v}^2$ (16)

where $m$ = The mass of the projectile (kg);

$g$ = Acceleration due to gravity (m/s²);

$v$ = The instantaneous velocity (m/s);

$c$ = Gravity coefficient (kg/m).

The term $k{v}^2$ represents the quadratic drag force, which becomes increasingly significant at higher velocities and is a more accurate representation of air resistance for many real-world objects compared to linear models.

If we consider the spherical projectile shape (baseball) thrown vertically upward with speed $v_0$ and is subjected to a quadratic drag with magnitude $f\left(v\right)=c{v}^2$ , and the velocity reaches its terminal ($v_{term}$ ), then the separation from the vertical equation of the motion becomes as in Equation (17).

$\frac{\text{d}{v}^{\prime }}{1-\frac{{v^{\prime }}^2}{v_{term}^2}}=gt$ (17)

In this separation form, we can integrate both sides directly (assuming $v_0=0$ the condition is met).

${\displaystyle \underset{0}{\overset{v}{\int }}\frac{\text{d}{v}^{\prime }}{1-\frac{{v^{\prime }}^2}{v_{term}^2}}=g{\displaystyle \underset{0}{\overset{t}{\int }}\text{d}{t}^{\prime }}}$ (18)

Consider the left-hand side of Equation (18):

${\displaystyle {\int }_0^v\frac{\text{d}{v}^{\prime }}{1-\frac{{v^{\prime }}^2}{v_{term}^2}}\Rightarrow v_{term}{\displaystyle {\int }_0^{\frac{v}{v_{term}}}\frac{\text{d}x}{1-x^2}}}$ $$ (19)

$\Rightarrow v_{term}{\text{arctanh}\left(x\right)|}_0^{\frac{v}{v_{term}}}\Rightarrow v_{term}\text{arctanh}\left(\frac{v}{v_{term}}\right)$ (20)

The simplification of (20) yields Equation (21) as the vertical velocity of the projectile motion under quadratic drag.

$v\left(t\right)=v_{term}\tanh \left(\frac{gt}{v_{term}}\right)$ (21)

To get the maximum vertical displacement of the projection under quadratic drag motion, we have to integrate the other form of equation $v=-g\left[1+{\left(\frac{v}{v_{term}}\right)}^2\right]$ , which yields Equation (22).

$y_{\max }\left(t\right)=\frac{v_{term}^2}{g}\ln \left(\frac{v_{term}-v_0}{v_{term}}\right)$ or $y_{\max }\left(t\right)=\frac{v_{term}^2}{g}\ln \left[\cosh \left(\frac{gt}{v_{term}}\right)\right]$ (22)

The Equations (21) and (22) characterize the vertical motion of a body subjected to gravitational acceleration in the presence of a quadratic drag force, where the resistive force is proportional to the square of the velocity. These models assume motion through a uniform medium with constant physical properties, neglecting all forces except gravity and quadratic drag on the form $F_{drag}=-c{v}^2$ . The velocity-time relation, given by Equation (21), describes the temporal evolution of velocity for a body released from rest, where $v_{term}=\sqrt{mg/c}$ denotes the terminal velocity attained as $t\to \infty$ . The associated displacement, expressed in Equation (22), represents the vertical position as a function of time for initial upward motion and symmetric ascent-descent scenarios, respectively. These analytical solutions are critical for accurately modelling projectile trajectories in real-world applications, particularly where air resistance cannot be neglected, such as in aerospace vehicle dynamics, sports ballistics, and environmental physics.

3. Results and Discussion

This section presents the computational results obtained by applying the fourth-order Runge-Kutta (RK4) method to simulate projectile motion under quadratic air resistance. The influence of varying projectile masses and initial velocities on horizontal and vertical kinematic parameters is systematically analysed. The results presented herein not only validate the underlying theoretical formulations but also annotate the non-trivial effects of drag-induced asymmetry in projectile trajectories, which are often underestimated in linear or drag-free models.

3.1. Horizontal Displacement

The horizontal displacement data expose a clear nonlinear relationship with both projectile mass and initial launch velocity. As shown in Table 2, increasing the initial speed from 10 m/s to 50 m/s results in substantial gains in horizontal range. However, the magnitude of these gains is disproportionately larger for heavier projectiles. Figure 1(c) shows that the horizontal displacement for the softball (0.189 kg) at 50 m/s is approximately 698.81 m, compared to just 209.92 m for the golf ball (0.046 kg) in Figure 1(a). This disparity accentuates the inverse scaling of the quadratic drag force with mass, as articulated by the ratio $\frac{c{v}^2}{m}$ , where lighter bodies experience more pronounced deceleration.

Table 2. The variations of initial velocities and masses with corresponding displacements (m).

Mass (kg)	Initial Velocity (m/s)
	50	40	30	20	10
0.046	209.9205	202.6346	193.2464	180.0276	157.4870
0.145	543.6866	520.8102	491.3668	449.9971	379.8297
0.198	698.8133	667.6405	627.5437	571.2683	476.0876

Figure 1. The horizontal displacement of different mass projectiles: (a) golf ball (0.046 kg); (b) baseball (0.145 kg); (c) softball (0.198 kg).

The simulations further demonstrate that at lower initial velocities, the displacement curves tend to flatten, indicative of a saturation behavior commonly observed in drag-dominated regimes. The diminished marginal gains in displacement at lower speeds reinforce the quadratic nature of the resistive force, which imposes exponentially increasing energy losses with velocity. This highlights the inadequacy of ideal projectile models in accurately predicting flight range in realistic aerodynamic environments, particularly for lightweight projectiles (see Figure 1).

3.2. Horizontal Velocity

Table 3 shows the horizontal velocity of each projectile mass after a flight duration of 40 seconds. Figure 2 shows the results of a steep decay in horizontal velocity for all cases, with the rate of decay inversely correlated with mass. At 50 m/s initial velocity, the golf ball’s horizontal speed drops to 0.80 m/s, while the baseball and softball maintain significantly higher terminal velocities of 2.45 m/s and 3.15 m/s, respectively. This is consistent with the analytical expectation that terminal velocity in the presence of quadratic drag is attained faster in lighter bodies due to their reduced inertial resistance to deceleration.

The decay trends observed are highly consistent with the analytical expression in Equation (11), where Equation (12) represents the characteristic decay time. The slower decay observed in heavier projectiles reflects their longer characteristic times, reaffirming the theoretical model and verifying the numerical accuracy of the simulation. Notably, the approach to terminal velocity is asymptotic, demonstrating the diminishing impact of acceleration over time due to drag dominance (Table 3).

Table 3. The final velocities (m/s) with mass variations at time reference 40 s.

Mass (kg)	Initial Velocity (m/s)
	50	40	30	20	10
0.046	0.8046	0.8014	0.7961	0.7857	0.7560
0.145	2.4514	2.4217	2.3738	2.2835	2.0495
0.198	3.1484	3.0996	3.0216	2.8767	2.5150

Figure 2. The horizontal final velocities of different mass projectiles: (a) golf ball (0.046 kg); (b) baseball (0.145 kg); (c) softball (0.198 kg).

3.3. Vertical Displacement

The vertical displacement of a projectile subjected to quadratic air drag was analyzed for varying projectile masses (0.046 kg, 0.145 kg, and 0.198 kg) over different time intervals. The results, presented in Figures 3(a)-(d) and summarised in Table 4, demonstrate the significant influence of both mass and time on the vertical displacement trajectory under drag effects.

As shown in Figure 3, the displacement-time curves disclose nonlinear behaviour due to the velocity-squared dependence of the drag force. For all cases, the projectile with the highest mass (0.198 kg) consistently attained the most significant vertical displacement, while the lightest projectile (0.046 kg) experienced the most pronounced reduction in displacement. This behaviour is expected, as the decelerating effect of drag is inversely related to mass, causing lighter objects to lose vertical momentum more rapidly.

It is observed that at 5 seconds, the vertical displacement ranges from 67.03 m (for 0.046 kg) to 97.94 m (for 0.198 kg). By 20 seconds, the displacement increases significantly, reaching 645.72 m for the 0.198 kg projectile. This progression, detailed in Table 4, further underscores the nonlinear growth of displacement and the mass-dependent nature of drag deceleration.

Moreover, the velocities shown in each subfigure highlight how drag reduces the effective vertical velocity more significantly for lighter projectiles. The declining slope gradient over time, especially for the 0.046 kg mass, confirms that energy dissipation due to drag becomes increasingly dominant as time progresses. The corresponding velocities at the end of each time interval are also indicated in the legends. The results illustrate that greater projectile mass leads to reduced sensitivity to drag force, resulting in higher vertical displacements over time.

Table 4. Vertical displacement (m) for varying projectile masses and time intervals under quadratic drag.

Mass (kg)	Time Variation
	5 s	10 s	15 s	20 s
0.046	67.03	156.46	246.03	335.60
0.145	92.16	246.79	405.60	564.62
0.198	97.94	274.77	459.93	645.72

Figure 3. Vertical displacement of a projectile under the influence of quadratic air drag for three different masses: 46 g (blue), 145 g (red), and 198 g (black). Subfigures (a) - (d) represent displacement profiles over increasing time intervals: (a) 0 - 5 s, (b) 0 - 10 s, (c) 0 - 15 s, and (d) 0 - 20 s.

These results validate the theoretical model of quadratic drag and suggest that mass is a critical parameter when predicting vertical motion in resistive media. The consistency between simulated trajectories and physical expectations indicates that the numerical model reliably captures drag-induced motion suppression.

3.4. Vertical Velocity

A consistent and central finding across both horizontal and vertical analyses is the pivotal role of mass in modulating drag effects. Heavier projectiles exhibit greater mechanical robustness against the retarding influence of air resistance, maintaining higher velocities and displacements throughout their trajectories. This is not merely a scaling effect; instead, it reveals a multiscale interaction between inertial and aerodynamic forces that profoundly alters trajectory characteristics. The differential acceleration profiles suggest that optimising projectile design for distance or stability in real-world applications must consider the nuanced interplay between mass, drag coefficient, and initial kinetic energy (Figure 4). Three masses 0.046 kg (blue), 0.145 kg (red), and 0.198 kg (black) are simulated, with corresponding terminal velocities of 17.9 m/s, 31.8 m/s, and 37.2 m/s, respectively. The curves clearly show how increased mass results in higher terminal velocity due to greater gravitational force overcoming the quadratic drag, consistent with the theoretical dependence $v_{term}\propto \sqrt{m}$ as in Equation (20).

Figure 4. Vertical velocity as a function of time for projectiles with different masses falling under gravity with quadratic air resistance.

Such mass-dependent behavior has implications beyond simple range estimation. In high Reynolds number flows, as simulated here, the projectile’s response to aerodynamic damping can inform predictive modeling in fields as diverse as aerospace re-entry, sports physics, and environmental dispersion of aerosols and debris.

3.5. Trajectory and Asymmetry Analysis

To directly visualize the combined effect of horizontal and vertical motion under quadratic drag and to assess the overall shape of the projectile path, we present the full 2D trajectories for the golf ball (0.046 kg), baseball (0.145 kg), and softball (0.198 kg) in Figure 5. All projectiles were launched from ground level with an initial velocity of 50 m/s at a launch angle of 45˚, which maximizes range in a vacuum and provides a clear baseline for comparison.

Figure 5. 2D trajectories for projectiles of different masses launched at 50 m/s and 45˚ under quadratic air drag.

The trajectories exhibit two key features predicted by theory but not visible in separate 1D analyses:

1) Trajectory Flattening: Compared to the ideal parabolic path (black dashed line, if added, or described as “vacuum trajectory”), all three trajectories are visibly flattened. The apex height is significantly reduced, and the range is shortened. This effect is most pronounced for the lightest projectile, the golf ball, whose peak height is less than half that of the vacuum trajectory. This demonstrates the dominant energy dissipation due to drag, particularly during the ascent phase when both gravity and drag act to decelerate the projectile.

2) Asymmetric Path: A crucial observation is the distinct asymmetry between the ascent and descent phases. The descent is steeper and more rapid than the ascent. This is because during ascent, the drag force acts downward, adding to gravity, resulting in a large negative acceleration. During descent, the drag force acts upward, opposing gravity, leading to a smaller net downward acceleration. This results in a trajectory that is skewed to the left of its apex compared to the symmetric parabola. The degree of asymmetry increases for lighter masses, as seen in the golf ball’s trajectory, which has a much sharper drop-off after its peak compared to the more gradual descent of the heavier softball.

This 2D visualization directly connects our component analyses to the overall projectile path. The reduced range observed in Table 4 and the lower maximum height noted in Section 3.3 are now seen as consequences of the distorted, asymmetric, and flattened trajectory caused by quadratic drag. The mass dependence is also evident: the heavier softball follows a path closest to the ideal parabola, while the lighter golf ball deviates most significantly, confirming that increased mass provides greater resistance to the distorting effects of aerodynamic drag

4. Conclusions

The fidelity of the numerical model $\theta = 45˚$ is emphasized by its consistency with theoretical predictions across multiple observable quantities, including displacement, velocity decomposition, terminal velocity, and trajectory asymmetry. The application of the RK4 method ensures numerical stability and accuracy, even in the presence of stiff differential equations arising from nonlinear drag terms. Notably, the computational approach retains physical interpretability, offering detailed insight into parameter sensitivities (e.g., mass, velocity, drag coefficient) that are otherwise obscured in purely analytical or overly simplified models.

The trajectory shapes obtained under quadratic drag differ qualitatively from ideal parabolic paths, exposing characteristic flattening and curvature distortions that are indicative of real-world aerodynamic interactions. This provides strong justification for incorporating realistic drag models into simulation frameworks, especially in high-speed or precision-dependent applications. Future work will include controlled experimental launches using high-speed video tracking to validate the predicted trajectories and terminal velocities, thereby establishing a critical empirical benchmark for the model’s predictive capabilities.

Acknowledgements

The authors acknowledge the support from Mbeya University of Science and Technology (MUST) and Abdulrahman Al-Sumait University (SUMAIT University).

Data Availability Statement

All data can be obtained from the corresponding authors upon request.

Conflicts of Interest

The authors declare no conflict of interest regarding the publication of this paper.

References

[1]	Linthorne, N.P. (2001) Optimum Release Angle in the Shot Put. Journal of Sports Sciences, 19, 359-372.[CrossRef] [PubMed]
[2]	Linthorne, N.P. and Everett, D.J. (2006) Soccer: Release Angle for Attaining Maximum Distance in the Soccer Throw-In. Sports Biomechanics, 5, 243-260.
[3]	Linthorne, N.P., Guzman, M.S. and Bridgett, L.A. (2005) Optimum Take-Off Angle in the Long Jump. Journal of Sports Sciences, 23, 703-712.[CrossRef] [PubMed]
[4]	El-Sayed, A., Nour, H., Raslan, W. and El-Shazly, E.S. (2015) A Study of Projectile Motion in a Quadratic Resistant Medium via Fractional Differential Transform Method. Applied Mathematical Modelling, 39, 2829-2835.[CrossRef]
[5]	Lubarda, M.V. and Lubarda, V.A. (2022) A Review of the Analysis of Wind-Influenced Projectile Motion in the Presence of Linear and Nonlinear Drag Force. Archive of Applied Mechanics, 92, 1997-2017.[CrossRef]
[6]	Bradshaw, J.L. (2023) Projectile Motion with Quadratic Drag. American Journal of Physics, 91, 258-263.[CrossRef]
[7]	Hackborn, W.W. (2016) On Motion in a Resisting Medium: A Historical Perspective. American Journal of Physics, 84, 127-134.[CrossRef]
[8]	Chudinov, P. (2014) Approximate Analytical Description of the Projectile Motion with a Quadratic Drag Force. Athens Journal of Sciences, 1, 97-106.[CrossRef]
[9]	Timmerman, P. and van der Weele, J.P. (1999) On the Rise and Fall of a Ball with Linear or Quadratic Drag. American Journal of Physics, 67, 538-546.[CrossRef]
[10]	Sadraey, M. and Müller, D. (2009) Drag Force and Drag Coefficient. In: Sadraey, M., Ed., Aircraft Performance: Analysis, VDM Verlag Dr. Müller, 46.
[11]	D’Auria, F., Camargo, C. and Mazzantini, O. (2012) The Best Estimate Plus Uncertainty (BEPU) Approach in Licensing of Current Nuclear Reactors. Nuclear Engineering and Design, 248, 317-328.[CrossRef]
[12]	Džiugys, A. and Peters, B. (2001) An Approach to Simulate the Motion of Spherical and Non-Spherical Fuel Particles in Combustion Chambers. Granular Matter, 3, 231-266.[CrossRef]
[13]	Thornton, S.T. and Marion, J.B. (2021) Classical Dynamics of Particles and Systems. Cengage Learning.
[14]	Wohlwill, E. (2001) The Discovery of the Parabolic Shape of the Projectile Trajectory. Science in Context, 14, 375-410.[CrossRef]
[15]	Ray, S. and Fröhlich, J. (2014) An Analytic Solution to the Equations of the Motion of a Point Mass with Quadratic Resistance and Generalizations. Archive of Applied Mechanics, 85, 395-414.[CrossRef]
[16]	Yabushita, K., Yamashita, M. and Tsuboi, K. (2007) An Analytic Solution of Projectile Motion with the Quadratic Resistance Law Using the Homotopy Analysis Method. Journal of Physics A: Mathematical and Theoretical, 40, 8403-8416.[CrossRef]
[17]	Nayfeh, A.H. (2024) Perturbation Methods. Wiley.
[18]	Patil, D. (2025) Computational Physics: Basic Concepts. Educohack Press.
[19]	Ramaswamy, J. (2025) Applications of Differential Equations. Educohack Press.
[20]	Wang, J. (2015) Computational Modeling and Visualization of Physical Systems with Python. Wiley.
[21]	Attaullah, Y., Alyobi, S., Al-Duais, F.S., et al. (2023) On the Comparative Performance of Fourth Order Runge-Kutta and the Galerkin-Petrov Time Discretization Methods for Solving Nonlinear Ordinary Differential Equations with Application to Some Mathematical Models in Epidemiology. AIMS Mathematics, 8, 3699-3729.[CrossRef]
[22]	Saltelli, A., Ratto, M., Andres, T., et al. (2008) Global Sensitivity Analysis: The Primer. Wiley.
[23]	Hereman, W. (2006) Symbolic Computation of Conservation Laws of Nonlinear Partial Differential Equations in Multi-Dimensions. International Journal of Quantum Chemistry, 106, 278-299.[CrossRef]
[24]	Said, A.A., Mshewa, M.M., Mwakipunda, G.C., Ngata, M.R. and Mohamed, E.A. (2023) Computational Solution to the Problems of Projectile Motion under Significant Linear Drag Effect. Open Journal of Applied Sciences, 13, 508-528.[CrossRef]
[25]	Sahu, J. (2008) Time-Accurate Numerical Prediction of Free-Flight Aerodynamics of a Finned Projectile. Journal of Spacecraft and Rockets, 45, 946-954.[CrossRef]
[26]	Cross, R. (1998) The Sweet Spot of a Baseball Bat. American Journal of Physics, 66, 772-779.[CrossRef]
[27]	Pujol, O. and Pérez, J.P. (2007) On a Simple Formulation of the Golf Ball Paradox. European Journal of Physics, 28, 379-384.[CrossRef]
[28]	Escamilla, R.F. and Andrews, J.R. (2009) Shoulder Muscle Recruitment Patterns and Related Biomechanics during Upper Extremity Sports. Sports Medicine, 39, 569-590.[CrossRef] [PubMed]
[29]	Kane, T.R. and Levinson, D.A. (1985) Dynamics, Theory and Applications. McGraw Hill.
[30]	Lubarda, M.V. and Lubarda, V.A. (2022) A Review of the Analysis of Wind-Influenced Projectile Motion in the Presence of Linear and Nonlinear Drag Force. Archive of Applied Mechanics, 92, 1997-2017.[CrossRef]
[31]	Camargo, J.D.R., Ruiz, J.P. and Roa, C.A.B. (2024) Dynamic Analysis of Rounded Projectiles: Software Solution Development. Ciencia y poder aéreo, 19, 46-57.
[32]	Glavelis, T., Ploskas, N. and Samaras, N. (2010) A Computational Evaluation of Some Free Mathematical Software for Scientific Computing. Journal of Computational Science, 1, 150-158.[CrossRef]
[33]	Sharma, N. (2010) A Comparative Study of Several Numerical Computational Packages. University of Maryland, Baltimore County.
[34]	Rowan, M. (2006) Probabilistic Analysis of Pulsed Signals in Radio Telescope Data Using Particle Filters. https://scholar.google.com/scholar?hl=en&as_sdt=0%2C5&q=Probabilistic+Analysis+of+Pulsed+Signals+in+Radio+Telescope+Data+using+Particle+Filters&btnG=
[35]	Yokogawa, T., Hannan, M.C. and Burgess, H.A. (2012) The Dorsal Raphe Modulates Sensory Responsiveness during Arousal in Zebrafish. The Journal of Neuroscience, 32, 15205-15215.[CrossRef] [PubMed]
[36]	Oberkampf, W.L. and Trucano, T.G. (2002) Verification and Validation in Computational Fluid Dynamics. Progress in Aerospace Sciences, 38, 209-272.[CrossRef]
[37]	Oberkampf, W.L., Trucano, T.G. and Hirsch, C. (2004) Verification, Validation, and Predictive Capability in Computational Engineering and Physics. Applied Mechanics Reviews, 57, 345-384.[CrossRef]
[38]	Romate, J.E. (1990) The Numerical Simulation of Nonlinear Gravity Waves. Engineering Analysis with Boundary Elements, 7, 156-166.[CrossRef]
[39]	Speičvcys, L., Jensen, C.S. and Kligys, A. (2003) Computational Data Modeling for Network-Constrained Moving Objects. Proceedings of the 11th ACM International Symposium on Advances in Geographic Information Systems, New Orleans, 7-8 November 2003, 118-125.[CrossRef]
[40]	Yilmaz, A., Javed, O. and Shah, M. (2006) Object Tracking: A Survey. ACM Computing Surveys, 38, 13.[CrossRef]
[41]	Farkas, Z. (2011) Considering Air Density in Wind Power Production. arXiv:1103.2198.
[42]	Richter, A. and Nikrityuk, P.A. (2012) Drag Forces and Heat Transfer Coefficients for Spherical, Cuboidal and Ellipsoidal Particles in Cross Flow at Sub-Critical Reynolds Numbers. International Journal of Heat and Mass Transfer, 55, 1343-1354.[CrossRef]
[43]	Castang, C., Laín, S., García, D. and Sommerfeld, M. (2022) Aerodynamic Coefficients of Irregular Non-Spherical Particles at Intermediate Reynolds Numbers. Powder Technology, 402, Article 117341.[CrossRef]
[44]	Phillips, R.L. and Cohen, C.B. (1959) Use of Drag Modulation to Reduce Deceleration Loads during Atmospheric Entry. ARS Journal, 29, 414-422.[CrossRef]
[45]	Petrushov, V. (1998) Improvement in Vehicle Aerodynamic Drag and Rolling Resistance Determination from Coast-Down Tests. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, 212, 369-380.[CrossRef]