Projectile motion: optimal launch angle for weak quadratic drag

In the absence of air resistance, a launch angle of 45° maximises range. When there is drag linear in speed, the equations of motion can be integrated analytically, and closed-form expressions in terms of the Lambert W-function can be obtained for the optimal launch angle; see Packel & Yuen (2004). However, a more realistic model of air resistance has drag proportional to the square of speed, for which the equations of motion are nonlinear and analytic solutions cannot be obtained.

In high school there was this rather horrible investigation where we had to (experimentally) compare the optimal launch angle of a golf ball and a ping-pong ball. At the time I believed that nothing much could be done in terms of modelling other than solving the equations of motion numerically, but back then I knew nothing of scaling and perturbation theory. Now armed with some basic knowledge of these very useful tools, I have been able to derive an asymptotic expansion for the optimal launch angle when air resistance is relatively weak.

Solution

Diagram for projectile motion with air resistance, which causes the trajectory to be asymmetric instead of parabolic.

Equations of motion

Suppose a projectile of mass mm is launched at speed uu and angle ϕ\phi from the ground, which has gravitational field strength gg, and let there be drag proportional to the square of the projectile's speed, with constant of proportionality bb. Let dots denote time derivatives. The drag force on the projectile has magnitude bx˙2b \norm{\dot{\vec{x}}}^2, and it acts in the direction opposite to the projectile's velocity, i.e. in the direction x˙/x˙-\dot{\vec{x}} / \norm{\dot{\vec{x}}}. Therefore the drag force is bx˙x˙-b \dot{\vec{x}} \norm{\dot{\vec{x}}}, and so in components the equations of motion are

mx¨=bx˙x˙2+y˙2,x˙(0)=ucosϕ,x(0)=0,my¨=mgby˙x˙2+y˙2,y˙(0)=usinϕ,y(0)=0.\begin{alignedat}{4} m \ddot{x} &= && - b \dot{x} \sqrt{\dot{x}^2 + \dot{y}^2}, \quad & \dot{x} (0) &= u \cos\phi, \quad & x (0) &= 0, \\ m \ddot{y} &= -m g && - b \dot{y} \sqrt{\dot{x}^2 + \dot{y}^2}, \quad & \dot{y} (0) &= u \sin\phi, \quad & y (0) &= 0. \end{alignedat}

Scaling

Since we shall be making a perturbation from the dragless b=0b = 0 case, it is appropriate to choose the length scale and time scale thereof. In the dragless case the mass mm is irrelevant, and the only parameters are the initial speed uu and gravitational acceleration gg, yielding the length scale

L=u2/gL = u^2 / g

and the time scale

T=u/g.T = u / g.

Therefore we put x^=x/L\hat{x} = x / L, y^=y/L\hat{y} = y / L, and t^=t/T\hat{t} = t / T to obtain scaled (dimensionless) variables x^\hat{x}, y^\hat{y}, and t^\hat{t}. Dropping the hats, the equations of motion become

x¨=Bx˙x˙2+y˙2,x˙(0)=cosϕ,x(0)=0,y¨=1By˙x˙2+y˙2,y˙(0)=sinϕ,y(0)=0,\begin{alignedat}{4} \ddot{x} &= && - B \dot{x} \sqrt{\dot{x}^2 + \dot{y}^2}, \quad & \dot{x} (0) &= \cos\phi, \quad & x (0) &= 0, \\ \ddot{y} &= -1 && - B \dot{y} \sqrt{\dot{x}^2 + \dot{y}^2}, \quad & \dot{y} (0) &= \sin\phi, \quad & y (0) &= 0, \end{alignedat}

where

B=bu2mgB = \frac{b u^2}{m g}

is the initial drag-to-weight ratio, the only dimensionless group in the problem. By definition, the projectile's terminal speed cc is given by

1=bc2mg.1 = \frac{b c^2}{m g}.

Dividing these, we see that

B=(uc)2.B = \roundbr{\frac{u}{c}}^2.

Now the optimal angle is dimensionless, so it can depend only on the sole dimensionless group BB. Thus, the optimal angle depends only on B=u/c\sqrt{B} = u / c, the ratio between the initial and terminal speeds. (I wish I knew this back in Year 12.)

Perturbed trajectory

By "weak drag" I mean that B1B \ll 1, i.e. u2c2u^2 \ll c^2. We make an asymptotic expansion in powers of BB about B=0B = 0:

x˙=x˙0+x˙1B+x˙2B2+ ⁣O(B3),y˙=y˙0+y˙1B+y˙2B2+ ⁣O(B3).\gdef\grey#1{\textcolor{595959}{#1}} \gdef\plusorder#1{\mathbin{\grey{+}} \grey{\order\roundbr{#1}}} \gdef\coeffbr#1{\Bigl(#1\Bigr)} \begin{aligned} \dot{x} &= \dot{x}_0 + \dot{x}_1 B + \dot{x}_2 B^2 \plusorder{B^3}, \\ \dot{y} &= \dot{y}_0 + \dot{y}_1 B + \dot{y}_2 B^2 \plusorder{B^3}. \end{aligned}

Substituting these into the equations of motion and equating coefficients, we obtain successive 2nd-order ordinary differential equations for x0x_0, y0y_0, x1x_1, y1y_1, etc. Thus x˙0\dot{x}_0, x0x_0, y˙0\dot{y}_0, y0y_0, x˙1\dot{x}_1, …, y2y_2 can be determined by straight integration, the details of which I leave to the manuscript (page 3 onwards).

You might be wondering why we have stopped at order B2B^2. Initially I thought that integration could be performed to arbitrary order (although the amount of algebra required grows very quickly). However I was wrong; it turns out that the following integrals appear at cubic order, which have no closed-form expression:

τlog(τ+α2+τ2)α2+τ2 ⁣dτ,log2(τ+α2+τ2)(α2+τ2)3/2 ⁣dτ.\begin{aligned} & \int \frac{ \tau \log \roundbr{\tau + \sqrt{\alpha^2 + \tau^2}} }{ \alpha^2 + \tau^2 } \td\tau, \\[\tallspace] & \int \frac{ \log^2 \roundbr{\tau + \sqrt{\alpha^2 + \tau^2}} }{ (\alpha^2 + \tau^2)^{3/2} } \td\tau. \end{aligned}

If you are able to evaluate either of the two integrals above, please contact me at yawnoc.outsell414@passmail.net. On the other hand, quartic terms, even if they could be found, would be of little practical use, since the expansion is asymptotic and as such probably does not converge.

Flight time

Having determined the trajectory, we then determine the flight time, given by y(t)=0y (t) = 0. To quadratic order, i.e. with

y=y0+y1B+y2B2+ ⁣O(B3),t=t0+t1B+t2B2+ ⁣O(B3),\begin{aligned} y &= y_0 + y_1 B + y_2 B^2 \plusorder{B^3}, \\ t &= t_0 + t_1 B + t_2 B^2 \plusorder{B^3}, \end{aligned}

this becomes

y0(t0)+(t1y˙0(t0)+y1(t0))B+(t2y˙0(t0)+12t12y¨0(t0)+t1y˙1(t0)+y2(t0))B2+ ⁣O(B3)=0.y_0 (t_0) + \coeffbr{t_1 \dot{y}_0 (t_0) + y_1 (t_0)} B + \coeffbr{ t_2 \dot{y}_0 (t_0) + \tfrac{1}{2} {t_1}^2 \ddot{y}_0 (t_0) + t_1 \dot{y}_1 (t_0) + y_2 (t_0) } B^2 \plusorder{B^3} = 0.

The unperturbed (dragless) flight time is given by the positive solution to y0(t0)=t0(sinϕt0/2)=0y_0 (t_0) = t_0 (\sin\phi - t_0 / 2) = 0, which is

t0=2sinϕ.t_0 = 2 \sin\phi.

From the linear and quadratic terms we then obtain t1t_1 and t2t_2 (see page 20 onwards of manuscript).

Range

Substituting the flight time into the horizontal component of the trajectory gives the range

R=x(t)=x0(t0)+(t1x˙0(t0)+x1(t0))B+(t2x˙0(t0)+t1x˙1(t0)+x2(t0))B2+ ⁣O(B3)=R0+R1B+R2B2+ ⁣O(B3),\begin{aligned} R &= x (t) \\ &= x_0 (t_0) + \coeffbr{t_1 \dot{x}_0 (t_0) + x_1 (t_0)} B + \coeffbr{ t_2 \dot{x}_0 (t_0) + t_1 \dot{x}_1 (t_0) + x_2 (t_0) } B^2 \plusorder{B^3} \\ &= R_0 + R_1 B + R_2 B^2 \plusorder{B^3}, \end{aligned}

where R0R_0, R1R_1, and R2R_2 are (horribly complicated) functions of ϕ\phi (see page 24 onwards of manuscript).

Optimal launch angle

Let primes denote ϕ\phi derivatives. Then the optimal launch angle is given by R(ϕ)=0R' (\phi) = 0, or, to quadratic order,

R0(ϕ0)+(ϕ1R0(ϕ0)+R1(ϕ0))B+(ϕ2R0(ϕ0)+12ϕ12R0(ϕ0)+ϕ1R1(ϕ0)+R2(ϕ0))B2+ ⁣O(B3)=0.R'_0 (\phi_0) + \coeffbr{\phi_1 R''_0 (\phi_0) + R'_1 (\phi_0)} B + \coeffbr{ \phi_2 R''_0 (\phi_0) + \tfrac{1}{2} {\phi_1}^2 R'''_0 (\phi_0) + \phi_1 R''_1 (\phi_0) + R'_2 (\phi_0) } B^2 \plusorder{B^3} = 0.

From the constant term we have R0(ϕ0)=2cos(2ϕ0)=0R'_0 (\phi_0) = 2 \cos (2 \phi_0) = 0, yielding the familiar

ϕ0=π4\phi_0 = \frac{\pi}{4}

in the absence of air resistance. From the linear and quadratic terms we obtain after much differentiation and algebra (pages 27–35 of manuscript) the constants ϕ1\phi_1 and ϕ2\phi_2, and hence the result:

Result

For small initial drag-to-weight ratios (or small initial-to-terminal kinetic energy ratios)

B=bu2mg=u2c21,B = \frac{b u^2}{m g} = \frac{u^2}{c^2} \ll 1,

the optimal launch angle has the asymptotic expansion

ϕ=π4+(3322+164log1+1/211/2)B+(13933840+815122log1+1/211/2+172048log21+1/211/2)B2+ ⁣O(B3),\begin{aligned} \phi = \frac{\pi}{4} & + \coeffbr{ - \tfrac{3}{32} \sqrt{2} + \tfrac{1}{64} \log \tfrac{1 + 1 / \sqrt{2}}{1 - 1 / \sqrt{2}} } B \\ & + \coeffbr{ - \tfrac{1393}{3840} + \tfrac{81}{512} \sqrt{2} \log \tfrac{1 + 1 / \sqrt{2}}{1 - 1 / \sqrt{2}} + \tfrac{17}{2048} \log^2 \tfrac{1 + 1 / \sqrt{2}}{1 - 1 / \sqrt{2}} } B^2 \\ & \plusorder{B^3}, \end{aligned}

or

ϕ=45°6.018°B+3.290°B2+ ⁣O(B3)=45°6.018°u2c2+3.290°u4c4+ ⁣O(u6c6).\begin{aligned} \phi &= 45 \degree - 6.018 \degree B + 3.290 \degree B^2 \plusorder{B^3} \\ &= 45 \degree - 6.018 \degree \frac{u^2}{c^2} + 3.290 \degree \frac{u^4}{c^4} \plusorder{\frac{u^6}{c^6}}. \end{aligned}

Numerical verification

In the following table we compare numerically computed optimal launch angles with those from the asymptotic result above:

BB Optimal ϕ\phi
Numerical Asymptotic Rel. error
    0 45.0° 45.0°  0
    0.1 44.4° 44.4°  0.005%
    0.2 43.9° 43.9°  0.04%
    0.3 43.4° 43.5°  0.1%
    0.4 43.0° 43.1°  0.3%
    0.5 42.6° 42.8°  0.5%
    0.6 42.2° 42.6°  0.8%
    0.7 41.9° 42.4°  1.2%
    0.8 41.6° 42.3°  1.7%
    0.9 41.3° 42.2°  2.4%
    1 41.0° 42.3°  3.2%
    2 38.8° 46.1° 19%
    3 37.3°
    4 36.2°
    5 35.3°
    6 34.6°
    7 34.0°
    8 33.5°
    9 33.0°
   10 32.6°
   15 31.1°
   20 30.1°
   50 27.1°
  100 25.2°
  200 23.5°
  500 22.9°
 1000 19.7°
10000 18.1°

Indeed the asymptotic expansion is very accurate for B<0.5B < 0.5 (or equivalently u/c<0.7u / c < 0.7).

Graph of optimal launch angle (phi) versus the dimensionless group B.

The true optimal launch angle is a decreasing function of BB. Thus, very crudely, the asymptotic expansion becomes useless when it stops decreasing, which occurs at

Bϕ12ϕ2=0.9.B \approx \frac{-\phi_1}{2 \phi_2} = 0.9.

Finally, here is a plot of the optimal launch angle in terms of B=u/c\sqrt{B} = u / c, the initial-to-terminal speed ratio:

Graph of optimal launch angle (phi) versus the initial-to-terminal speed ratio (u divided by c).

See also

Cite this page

Conway (2023). Projectile motion: optimal launch angle for weak quadratic drag. <https://yawnoc.github.io/math/projectile-weak-drag> Accessed 2025-04-05.