Analytic Approximations of Projectile Motion with Quadratic Air Resistance

We study projectile motion with air resistance quadratic in speed. We consider three regimes of approximation: low-angle trajectory where the horizontal velocity, u, is assumed to be much larger than the vertical velocity w; high-angle trajectory where w u  ; and split-angle trajectory where w u  . Closed form solutions for the range in the first regime are obtained in terms of the Lambert W function. The approximation is simple and accurate for low angle ballistics problems when compared to measured data. In addition, we find a surprising behavior that the range in this approximation is symmetric about / 4  , although the trajectories are asymmetric. We also give simple and practical formulas for accurate evaluations of the Lambert W function.


Introduction
In a previous paper on projectile motion with air resistance linear in speed, we presented closed form solutions for the range in terms of the Lambert W function [1].Amid the growing list of problems that benefited from using the W function [2][3][4][5][6][7], one question naturally arises as to whether a similar approach exists if air resistance is quadratic in speed, the more realistic case in practice.
We have studied this problem and found that solutions exist for low-angle trajectories using the W function.In this paper, we report our findings, starting with an overview of the regimes of approximation in Section 2. We focus on the low-angle regime in Section 3 and discuss the dynamics, which leads to a remarkable property that the range is symmetric about / 4  , even in the presence of air resistance.For completeness, high-angle and splitangle regimes are briefly discussed in Section 4 and 5, respectively, followed by a comparison with observed data and discussions in Section 6.To make the W function easily accessible, in the Appendix we give simple and practical formulas for the accurate evaluation of this function.

Regimes of Approximation
We assume the net force F  , including air resistance and gravity on the projectile of mass, m, to be where b is the drag coefficient with the dimensions m -1 .The components of the velocity are v ui wj      , with / u dx dt  and / w dy dt  , where x and y are the usual horizontal and vertical positions of the projectile.The initial position of the projectile is at the origin.Throughout the paper, we use the following notations for frequently occurring terms: u 0 , w 0 = initial horizontal and vertical velocities; R 0 = 0 0 2u w g = range with no air resistance; To our knowledge, closed form solutions to (1) are known only for special initial conditions, not for arbitrary initial conditions [8].The difficulty is with the coupling of u and w in the speed which makes the problem inseparable.Further approximations are necessary to solve (1) analytically.
We consider three approximations aimed at linearizing the speed according to the relative magnitudes of u and w: low-angle trajectory (LAT) approximation, u w  ; highangle trajectory (HAT) approximation, w u  ; and split-angle trajectory (SAT) approximation, u w  , as, , Each of the approximations will be discussed below, with the emphasis on LAT that is important in practice and in ballistics.

The Trajectory
We assume in this case the horizontal velocity is on average much greater than the vertical velocity, u w  .This will be the case if the firing angle is small, a case also discussed by Parker [9].Here the speed may be approximated as v u  according to (3), so that the equations of motion in the LAT approximation are One can solve for u first, after which w, x, and y can be obtained.Hereafter, we will omit non-essential intermediate steps.The solutions may be verified by substitution into the equations of motion.The solutions are .For comparison, we also show the trajectory from ideal projectile motion with no air resistance and the trajectory from the solutions with the full 2  v air resistance, Equation 1.We will simply refer to the former as the ideal motion, and the latter as the full solution.
The full solution is carried out by numerically integrating the equations of motion with the full 2 v resistance (1), using the Runge-Kutta method.(The two curves labeled HAT and SAT in Figure 1 are discussed in Sections 4 and 5.) The agreement between the LAT approximation and the full solution is good at 20 o (nearly indistinguishable in Figure 1(a)), and it becomes worse for larger angles.This is as expected since the assumption was that LAT is valid only at small angles.The range is much reduced compared to the ideal motion.Air resistance introduces in the trajectories a well-known backward-forward asymmetry.The ascending part of the trajectory is shallower and longer, and the descending part is steeper and shorter.

The Range and Height in Comparison with Measurements
To find the range of the projectile we eliminate t from (6) to obtain The range, R, is the value of x when 0. y  It satisfies This is a transcendental equation that until now had been customarily solved numerically or graphically [9].But as we reported [1], equations of this type can also be solved analytically, with the results written in closed form in terms of the Lambert W function.
The Lambert W function [10] is defined, for a given value z, as the (inverse) function satisfying Our discussions below refer closely to the properties of this relatively "new" function.For readers unfamiliar with this function, we give a brief review in the Appendix.There the reader will also find some practical formulas for evaluating W.
To solve for the range R in terms of W, (8) needs to be rearranged in the form of ( 9) such that the multiplicative prefactor to the exponential,   W z , is the same as the exponent.This can be achieved by following the steps in [1].The result is (10) with defined in (2).We identify from ( 9) and ( 10) that This is the closed form expression for the range, R, in the LAT approximation.
The height, H, can be obtained by maximizing y in (7), and is given by We compare in Table 1 measured range and height data with calculations from Equations 11 and 12 for firing angles less than one degree.(The details of the calculations are given in the next subsection.) Table 1 shows that the measured data [11] and the calculations agree well once the firing angle is above about 10 minutes of arc.The discrepancy between measurement and theory is due to the uncertainty in the value b (see Table 1 caption), and not the approximation itself.Because the largest angle is still less than 1, higher order corrections are negligible.The relative error for the height is usually much larger than the relative error for the range.This is probably due to the difficulty in accurately measuring the relatively small height in this case.
The large errors in height and in range at the two smallest angles need not cause concern.It is due to the combination of exceeding difficulty in determining the small angle and the small height at the same time.Note that the diameter of the projectile and the height are of the same order of magnitude here.Overall, this example shows that if a reasonable b can be obtained, the LAT approximation should work well for low angle ballistics problems.

The Symmetry of Range in Firing Angle
The analytic Formula 11 enables us to immediately draw several surprising conclusions on the general properties of the range, R. Since z is a function of  which depends on the product we conclude that i) the range R is symmetric about / 4  , that is to say, two firing angles 1  and 2  will lead to  Given initial conditions 0 u and 0 w , the range, R, can be computed from (11).Since z is negative and W(z) is multi-valued (see Appendix) for 0 z ＜ , we have a choice of the branches 0 W or 1 W  .Comparing z in ( 11) and ( 9) we identify one trivial solution, namely the primary branch,

(see the LAT curves in
This solution, although mathematically correct, is unphysical because it gives a zero range when substituted into (11).The physical choice must be 1 W  .We note that for linear resistance [1] the choices were the opposite, where 0 W was the physical solution and 1 W  the unphysical one.
The range, R, calculated from (11) using  The LAT range is in good agreement with the full solution for low firing angles as expected, up to around / 6  .Compared to the ideal case, the maximum ranges in the LAT and the full solutions are substantially reduced, by about 40% for this particular set of parameters.We note that the LAT approximation produces asymmetric trajectories (Figure 1) but symmetric ranges.The physical reason can be traced to two factors influencing the range: the time of flight and the average horizontal velocity, as discussed below.

The Time of Flight and the Average Velocity
The average horizontal velocity, u , and the time of flight, T, are related by .

R u T
The time of flight T can be obtained from ( 6) by setting x R  , at t T  .Together with (11), this gives The average horizontal velocity can be expressed from ( 14) and ( 15) as Quantitatively, as the firing angle increases, the time of flight T increases, but the average horizontal velocity u decreases.However, before / 4  , the increase in T is more than the decrease in u so that the range as governed by ( 14), increases.After / 4  , however, the opposite happens so that the range decreases.The symmetric range is a result of the balance between T and u .

The Range for Small and Large Air Resistance
In the limit of small air resistance, 0 b  , the dimensionless parameter  will be small, 0 , and Using the properties of W(z) and after some algebra (we leave the details as an exercise to the interested reader), the first order correction to the range, R, in (9) is where 0 R is the range of ideal projectile motion, 0 b  .For large air resistance, 1 b  and 1 according to (11).Using the as-ymptotic expressions for (see Appendix), we have from (11)   It is interesting to compare the scaling behavior with our earlier study [1] for linear resistance, where we found the range to scale as 1/ b , the inverse of the resistance.For quadratic resistance, (18) indicates This shows that the logarithm term is characteristic of the quadratic resistance.

High-Angle Trajectory (HAT) Approximation
When the firing angle is large (close to / 2  ), we expect that, on average, the vertical velocity w will be much larger than the horizontal velocity, w u  .The speed is approximated as v w  from (3).The equations of motion in the HAT approximation are The solutions are broken into two parts because of |w|.In the ascending part of the trajectory, the solutions are With the values of u, w, x, y in (20,21) at the top as the initial condition for the descending trajectory, the solutions for descent are The time  starts from zero (at the top) in (22,23).The trajectories in the HAT approximation are shown in Figure 1 The best agreement with the full solution is seen at the highest angle 70 o , consistent with the underlying assumptions.We note that the agreement is considerably worse descending than ascending (Figure 1, 70 o (c)).The reason is that near the top, 0 w  , and the validity of the HAT approximation breaks down, causing the large discrepancy while falling back down.By contrast, the LAT approximation (Figure 1, 20 o (a)) is valid globally as long as the firing angle is small, giving a much better agreement on both parts of the trajectory.

Split-Angle Trajectory (SAT) Approximation
In Sections 3 and 4 we discussed low and high angle trajectories.To be complete, we consider in this section the split angle / 4 Note that upon replacing 2b b  in (24), / du dt is the same as that in (4) of LAT, and / dw dt is the same as that in (19) of HAT.The solutions for u, x will be the same as those in (5,6), and the solutions for w, y will be the same as for w, y in (21,23), so they will not be repeated here.
Similarly, the trajectories can be computed as before (with b replaced by 2b , of course).They are also shown in Figure 1 at the same angles with the same parameters.Here, we see the best agreement with the full solution at 45 o as it should.But, unlike the LAT or HAT curves, where after certain point in time (just before reaching the top) the differences keep increasing, the SAT curve crosses the full solution during the course of motion.This is due to the balance of the horizontal and vertical resistance forces.
Because of this balance, the SAT behaviors are interestingly different at low versus high angles.At 20 o (Figure 1(a)), the SAT curve is "squeezed" horizontally in comparison with the full solution, resulting in a shorter range and a higher height.This is because for low angles where u w ＞ , the horizontal resistance is over-estimated in the equations of motion (24).Conversely, at 70 o (Figure 1(c)), the SAT curve is compressed vertically, causing a lower height but a longer range.The reason is similarly due to the over-estimation of the resistance in the vertical direction.As a result, curve crossing occurs.

Conclusions
In summary, we have presented a detailed discussion of projectile motion with quadratic air resistance in three approximations.Our focus was on the low-angle trajec-tory approximation where we found closed form solutions for the range and the time of flight in terms of the secondary branch of the Lambert W function, 1 W  .The approximation is simple and accurate for low angle ballistics problems.
Various analytic properties were readily analyzed with these solutions.Together with projectile motion with linear air resistance [1,5], the example studied here serves two educational purposes: i) It is possible to introduce the use of special functions in physics at early undergraduate levels in a familiar, more realistic problem; and ii) It represents a good complimentary case where the physical solution required the secondary branch, 1 W  , rather than the principal branch 0 W as in the case of linear air resistance.
One interesting and closely-related question remains open, i.e., whether both branches of W might be required simultaneously in a physical solution, say when the air resistance contains both linear and quadratic terms, , under some forms of approximation, presumably.0 dR d  , from [1], yielding .
This implies either ,

Appendix 1. The Lambert W function
In this appendix, we briefly summarize some relevant elements of the Lambert W function.The reader can find an extensive review in [10].In lieu of the growing interest in this function in the physics community [2][3][4][5][6][7], we also present several practical formulas for the evaluation of the W function.
The Lambert W function as defined by ( 9) is generally complex-valued.For our purpose we are interested in the real-valued W(x), namely the principal branch 0 W , and the secondary branch 1 W  .The two branches are shown in Figure 3.
The behaviors of 0 W and 1 W  for small and large x are

Evaluation of W
During the course of our study, we were unable to find in the literature simple and practical algebraic expressions for evaluating W over the whole range, and for embedding it in our code.We therefore devised several readily usable approximate formulas given here for this purpose.It should give the interested reader a good starting point in using this function.
W is not (yet) available on a calculator like some elementary functions.It is debatable whether it should be elevated to the status of an elementary function.(See [13] for an enchanting account.)However, one can easily implement it on a programmable calculator with the formulas given here.The accuracy is at least 8 digits, in fact it is usually much better.

1) W in the Regular Regions
We first give the formulas optimized in the regular regions.Owing to space limitations, we will explain the methods used elsewhere [14].The general form of the expressions is where C is a constant and   The expansion variable, r, is related to the independent variable, x.We give in Table 2 the constant C, the variable r, and the coefficients i a and i b .To utilize the table, locate which function and region to use, then evaluate (A.2) with the corresponding C, r, i a , and i b .

2) W in the Asymptotic Regions
Both 0 W and 1 W  have the same form in the asymptotic regions

W
x q q q r a r a r a r a r p q x x p q r x p q a a q a q q a q q q   W  over all regions of x can be found by using Table 2 plus (A.4) for fixed precision (8 digits or better), or (A.5) and (A.6) for arbitrary precision.

6 )
Equation 6 gives the trajectory in the LAT approximation.The trajectories computed from (6) are plotted in Fig- ure 1 for three angles: 20 o , 45 o , and 70 o .The initial speed is 9.8 m/s for all angles.The drag coefficient is

Figure 1 .
Figure 1.The trajectories for an initial speed of 0 9.8m/s v  and drag coefficient 1 0.1m b   at three firing angles, 20 o (a), 45 o (b), and 70 o (c).The labeled curves are (see text): Full (thick solid line) -the numerical solution with the full 2 v resistance; LAT (solid line); HAT (dashed line, see Section 4); SAT (dash-dotted line, see Section 5), Ideal (dotted line) -motion without air resistance

Table 1 .
Measured and calculated results for a projectile of mass m = 9.7 g, diameter d = 0.76 cm, and muzzle velocity 823 m/s.The measured data are taken from [11].Results are calculated from the LAT approximation (11,12).The drag coefficient b = 1.05x10 -3 m -1 is determined from 2 / b C d m   , where C = 0.15 is the recommended value [11], and  = 1.2 kg/m

Figure 2 4 
as a function of the firing angle.It clearly demonstrates that R is symmetric about, and maximum at, / 4  .Also shown in Figure2, for com- parison, are the ranges for ideal projectile motion with no air resistance and with the full2  v resistance.Unlike the LAT case, the range of the full solution in Figure2is not symmetric about / 4  .It reaches maximum at an angle below / , a fact that is well known.

Figure 2 . 4 
Figure 2. The range of projectile motion as a function of the firing angle 0  .The initial speed is 0 v = 9.8 m/s and the drag coefficient is b = 0.1 m -1 .The LAT approximation (solid line) is maximum at and symmetric about / 4  , as is the ideal motion (dash-dotted line).The range with the full 2 v resistance (dashed line) peaks before / 4 

Figure 3 .and the secondary branch 1 W
Figure 3.The Lambert W function for real argument x.The two real branches are the principal branch 0 W (solid line), and the secondary branch

4 ) 1 W
This expression (A.4) should be used for 0 W in the rethat it is the subtle difference in p, namely,  that automatically selects the correct branch.

3 ) 5 )
W with a Little ProgrammingIf the reader wishes to calculate W with arbitrary precision, one can use Newton's rule[15] which, for a given x, is the root, w, to  The rootfinding process is as follows: Starting with an initial guess, 1 w , the successive iterations, n w , approach rapidly to the true value as It only remains to determine the initial guess, 1 w .One way to choose 1 w is

3 is the air density
Remarkably, these properties i) and ii) with air resistance in the LAT approximation are exactly the same as for ideal projectile motion without air resistance.

Table 2 . The approximate formulas for the evaluation of W with (A.2) and (A.3). Spaces are inserted in the coefficients for readability. For optimal precision, all the digits should be used Function
[14]an be shown[14]that (A.5) with the seed from (A.6) always converges to the correct branch 0