Projectile Motion with Resistance: Symbolic Solutions, the Lambert W Function, and the Inverse Range Problem

Ed Packel (packel@lfc.edu) and Dave Yuen (yuen@lfc.edu), Lake Forest College

Introduction

Mathematical investigations of projectile motion have a rich and vital history, going back almost 500 years. Besides the obvious application to ballistics, there is a much more noteworthy connection (for our more pacific purposes) to developments and colorful personalities in mathematics and physics. Many of these ideas are presented in a compelling paper by Groetsch [4], who traces a rich history from Tartaglia to Galileo and then employs a tour de force of undergraduate analysis to answer and expand some classical questions to the case of projectile motion in a resistive medium.

The purpose of this paper is to indicate how some of Groetsch's ingenious analysis can be circumvented with the help of computer algebra and a relatively new symbolic object, the Lambert W function, that increasingly seems destined for fame and immortality (see FOCUS newsletter [1]; Corless et al [3]). We will use these modern tools to simplify the derivation of some of Groetsch's results while extending some of them. In particular, we find a symbolic solution for the range as a function of the elevation angle in the presence of a linear resistance and we also give a partial solution to the inverse problem of finding the elevation angle(s) that give rise to given range values by obtaining a closed form for the angle generating the maximum range in terms of the initial velocity and the resistance constant.

Two interesting side issues arise in the process of our development. One of them evolves from the need to investigate certain limits involving the Lambert W function. To do this we develop a general theorem, which may be of interest in its own right, about inverse functions arising from real-analytic functions. A second issue relates to the newly emerging area of Experimental Mathematics (Borwein and Corless [2]) and raises practical and philosophical questions about use of symbolic computation to "discover" new results and the extent to which such computation can be viewed as an accepted form of proof.

Finding the Range With No Resistance; Struggling With Resistance

Deriving a formula in the absence of resistance for the range R as a function of the angle of elevation θ is a simple exercise in calculus. We race through it here to warm up for the more resistant case that follows. Consider a projectile that starts at the origin and is shot at an angle θ with an initial velocity v as pictured in Figure 1.

Figure 1.  Projectile motion with elevation angle θ and initial velocity v

Working with horizontal and vertical accelerations, we get the simple uncoupled pair differential equations.

		x'' = 0
y'' = -g.
		
Integrating each of these equations twice, using our initial conditions and noting that x' (0) = v cos(θ) and y' (0) = v sin (θ), with initial values for x and y both 0, we get

		x(t) = v cos(θ)*t
y(t) = ((-1)/2) g (t^2) + 
          v sin(θ)*t.
		
There are no surprises here -- it is easily shown that these parametric equations result in a parabolic path for the projectile. 

To get an equation for the range, we set the height y(t) to 0, compute the nonzero solution for t, and substitute it in the x(t) equation. This results in

		\(R = \(\(\(2 
                       v\^2\)\/g\) \(sin(θ)\)*\(cos(θ)\) = \(v\^2\/g\) \(sin(2  θ)\)\)\). 
		
As the plot of R vs. θ (we set v = 50and g = 32.2) in Figure 2 indicates, the maximal range occurs when θ = (π/4) and submaximal ranges occur for a pair of θ values equally spaced around (π/4). These results are independent of v and g and can be verified analytically using the symmetry of sin(2  θ) about θ = (π/4).

Figure 2.  A plot of range versus elevation angle

An attempt to extend these results to a projectile with resistance proportional (via constant k) to the velocity has a propitious beginning. The differential equations now become:

		x'' = -k x'
y'' = -g - k y'.
		
One integration on each equation (set u = x'and coast; set w = y', separate, and solve) leads (with initial values x' (0) = v cos(θ) and y' (0) = v sin(θ))

		x' = 
          v cos(θ) (e^(-k t ))
y'=\(\(\(\(-g\) + \((g + 
                          k\ v\ \(sin(θ)\))\)\ e\^\(\(-k\)\ t\)\)\/k\)\(.\)\)
		
		
		
Integrating again (with initial values x(0) = 0 and y(0) = 0) leads to

		\(x(
              t) = \(v\ \(cos(θ)\) \((1 - e\^\(\(-k\)\ t\))\)\)\/k\)
\(y(
              t) = \(\(\(\(-k\)\ t\ g + g + 
                    k\ v\ \(sin(θ)\) - e\^\(\(-k\)\ t\)\ \((g + k\ v\ \(sin(θ)\))\)\)\/k\^2\)\(.\)\)\)

Recall that in the "no resistance" case we set y(t) = 0, found a nonzero solution for t, and evaluated x(t) at this time t to express the horizontal range as a function of θ,  g,  and  v. 

Looking at the form of y(t) that results when resistance is present, we see that finding a nonzero root may be a daunting, if not impossible, task. Indeed, the presence of t and (e^(-k t))in an expression does not bode well for isolating t. The stage is now set for a dramatic rescue, so let us introduce the new function which will come to our aid.

The Lambert W Function and a Symbolic Solution for the Range

The Lambert W function can be defined as an inverse of the function T(w) = w (e^w) (see Figure 3).

Figure 3. The graph of w (e^w)

A look at the graph of T indicates that there are two such inverses (Figure 4), denoted by \(TraditionalForm\`W : \([\(-\(1\/e\)\), ∞)\) -> \ \([\(-1\), ∞)\)\)  and (Subscript[W, -1]) : [-(1/e), 
            0) ->  (-∞, -1].

Figure 4.  Two inverse functions for T(w) = w(e^w)

The evolution of W began with ideas proposed by J.H.Lambert in 1758 and refined by Euler over the subsequent two decades. An influential paper by Corless et al [3] presents a variety of applications, pure and applied, in which W plays a valuable role. As we will see, the range problem for a projectile with linear resistance is another such application.

By its inverse function definition, it follows directly that w = W(x) is a solution to the equation x = w (e^w). As a direct consequence, we have the following results:

	(1)   x = W(x) (e^(W(x)))  for all \(TraditionalForm\`x ∈ \([\(-\(1\/e\)\), ∞)\)\)
	
	(2)  W' (x) = ((W(x))/(x(1 + W(x))))   (follows neatly from the formula for the derivative of an inverse function and is a nice calculus exercise).
	
	(3)   A solution for t in the equation a t + b + 
          c (e^(  d t)) = 0 (with a d ≠ 0) is given by t=-(b/a)-(1/d)\(\(W(\(c\ d\ e\^\(-\(\(b\ d\)\/a\)\)\)\/a)\)\(.\)\)
	
To derive (3), we massage the first equation as follows:

	(a t  + 
              b) (e^(-d t)) = -c 
	
	\(TraditionalForm\`\((\(-d\)\ t - \(d\ b\)\/a)\) e\^\(\(-d\)\ t\) = \(c\ d\)\/a\)
	
	\(TraditionalForm\`\((\(-d\)\ t - \(d\ b\)\/a)\) e\^\(\(-d\)\ t - \(d\ b\)\/a\) = \(\(c\ d\)\/a\) e\^\(-\(\(d\ b\)\/a\)\)\).
	
Setting u = -d t - ((d b)/a), we then have \(TraditionalForm\`u\ e\^u = \(\(c\ d\)\/a\) e\^\(-\(\(d\ b\)\/a\)\)\)from which it follows that u = 
            W(((c d)/a) (e^(-((d b)/a)))). Substituting back for u and solving for t gives the desired result. 

Returning to the height equation y(t) = 0 in the form of the equation \(TraditionalForm\`\(-k\)\ t\  + 
          1 + \(\(\(k\)\(\ \)\(v\)\(\ \)\)\/g\) \(sin(θ)\) - e\^\(\(-k\)\ t\)\ \((1 + \(k\ v\)\/g\ \(sin(θ)\))\) = 
        0\), we substitute u = -1 - ((k v)/g) sin(θ) to get the simpler equation -k t - 
          u + (e^( -k t)) u = 0. Applying result (3) now gives the solution

	t = (1/k) (-u + 
                W(u (e^( u)))). 	
	
Before we find the range function, we deal with an interesting side issue. Since W is the inverse to T(w) = w (e^w), it would seem that applying W to both sides of this equation should give the identity w = W(w (e^w)). This is indeed a valid identity, but only for w ≥ -1since we restricted the domain of T to obtain its inverse W. But a look at the form of u = -1 - ((k v)/g) sin(θ) reveals that this quantity will never exceed -1 (constants are positive and θ ∈ [0, (π/2)]). As a result, our expression for t does not simplify, though if it did it would simply yield the obvious t = 0 solution, which is no help. It should be noted that \(TraditionalForm\`w\ e\^w ≥ \(-\(1\/e\)\)\) for all real w and is always in the domain of W.

Substituting our t value into x(t) gives us the desired range function, \(R(θ) = \(1\/k\) 
              v\ \(cos(θ)\) \((1 - \ e\^\(\(\ \)\(u - \ W(u\ e\^u)\)\))\)\).  From identity (1) above we may replace (e^( -W(
                  u (e^u)))) by \(\(W(
                u\ e\^\(\(\ \)\(u\)\))\)\/\(u\ e\^u\)\) since \(TraditionalForm\`u\ e\^\(\(\ \)\(u\)\) ≥ \(-\(1\/e\)\)\) is always true. This gives the somewhat simpler formula 

	\(R(θ) = \(1\/k\) 
              v\ \(cos(θ)\) \((1 - \(W(u\ e\^\(\(\ \)\(u\)\))\)\/u)\)\). 
	
Substituting for u, we get the range formula in all its elemental glory:

	\(R(θ) = \(1\/k\) 
              v\ \(cos(θ)\) \((1 - \(W(\((\(-1\) - \(\(k\ v\)\/g\) \(sin(θ)\))\)\ e\^\(\(\ \)\(\(-1\) - \(\(k\ v\)\/g\) \(sin(θ)\)\)\))\)\/\(\(-1\) - \(\(k\ v\)\/g\) \(sin(θ)\)\))\)\).

Computer Algebra Interlude: Honesty compels us to admit at this point that the idea for using Lambert W to find a closed form solution was really Mathematica's and not ours. Here is a sequence of commands that gave us a first closed form for the range.

x[t_] := 
        v*Cos[θ] (E^(-k t));
y[
          t_] := ((-k t g  + g + k v Sin[θ] - (E^(-k t)) (g + k v Sin[θ]))/(k^2));
tmax = Solve[y[t] == 0, t]

Simplify[x[t] /. tmax[[1]]]

Note that "ProductLog" is the Mathematica notation for Lambert W and that some algebra transforms the last output to the first form of R(θ) obtained above. It is somewhat amusing that Mathematica managed to lose the trivial "t = 0" root for y(t).

As a reality check on our range function, we plot in Figure 5 some sample graphs for varying resistance values k (keeping v = 50 and g = 32.2).

Figure 5.  Range versus angle of elevation for different k values

Note that, as one might expect, the maximal range angle decreases from (π/4) as resistance increases. Also, the symmetry about the maximal range angle for smaller ranges that was present in the no-resistance case disappears dramatically.

A Limit Theorem for a Class of Inverse Functions and Its Application

The theorem that follows was initially motivated by efforts to show that the range formula for R(θ) given above reduces to the classical "no resistance" result as the resistance k goes to zero. We will apply it to the Lambert W function for this purpose. Since the theorem holds for a much wider class of functions, we state and prove it in a more general context.

Theorem 1.   Let f(x) be a nonconstant real-analytic function in a neighborhood of (Subscript[x, 0]).  Suppose f(x) has a local extremum at (Subscript[x, 0]).  Then there is a neighborhood of (Subscript[x, 0]) such that f(x) is strictly monotonic on  [(Subscript[x, 0]), (Subscript[x, 0]) + h] and on [(Subscript[x, 0]) - h, (Subscript[x, 0])].  Let  g(x) be the inverse function of f(x) on [(Subscript[x, 0]), (Subscript[x, 0]) + h].  Then y(x) = g(f(x)) is defined on  [(Subscript[x, 0]) - k, (Subscript[x, 0])] for some   k > 0and \(lim\+\(x -> \(x\_0\^-\)\)\(y(x) - x\_0\)\/\(x - x\_0\)\  = \ \(-1. \)\) Note that this is the left-side derivative of y(x) because y((Subscript[x, 0])) = (Subscript[x, 0]).

Proof.  We have f(
          x) = (∑_(i = 0)^∞)(( (Subscript[a, i])(x - (Subscript[x, 0])))^i) in a neighborhood of (Subscript[x, 0]).  Let N be the least index N > 0 such that (Subscript[a, N]) ≠ 0 (such indices exist because f(x) is nonconstant).  So f(x) = 
        (Subscript[a, 0]) + (∑_(i = N)^∞)(( (Subscript[a, i])(x - (Subscript[x, 0])))^i).  By noting that f(x) = (Subscript[a, 0]) + (((x - (Subscript[x, 0])))^N) F(x), where \(TraditionalForm\`F(
          x) = ∑\+\(i = N\)\%∞\(\( a\_i\)(x - x\_0)\)\^\(i - N\)\), we see that F(x) has the same radius of convergence as f(x) and as a consequence must be continuous at (Subscript[x, 0]).  Since F((Subscript[x, 0])) = (Subscript[a, N]) ≠ 0, F(x) has constant sign in a neighborhood of (Subscript[x, 0]) and it is straightforward to see that N must be even in order for f(x) to have a local extremum at (Subscript[x, 0]).  Since the derivative f' (x) is real-analytic and not  identically zero, it is nonzero in some deleted neighborhood of (Subscript[x, 0]).  This implies that f(x) is strictly monotonic on the left and right parts of this deleted neighborhood separately. Using the notation y = y(x), we have, for any x ∈ [(Subscript[x, 0]) - k, (Subscript[x, 0])], f(y) = f(g(f(x)))=f(x)  since f(x) ∈ [(Subscript[x, 0]), (Subscript[x, 0]) + h] and f∘gis the identity function on this domain. Hence (Subscript[a, 0]) + (((y - (Subscript[x, 0])))^N) F(y) = 
        (Subscript[a, 0]) + (((x - (Subscript[x, 0])))^N) F(x) and \(TraditionalForm\`\((y - x\_0)\)\^N\/\((x - x\_0)\)\^N = \(F(x)\)\/\(F(y)\)\), Taking Nth roots of both sides and noting that ((y - (Subscript[x, 0]))/(x - (Subscript[x, 0]))) is negative, we have \(TraditionalForm\`\(y - x\_0\)\/\(x - x\_0\) = \(-\@\(\(F(x)\)\/\(F(y)\)\)\%N\)\). The desired result, \(lim\+\(x -> \(x\_0\^-\)\)\(y(x) - x\_0\)\/\(x - x\_0\)\  = \ \(-1\)\), follows directly by bringing the limit inside the radical on the right hand side.

By applying this theorem to the function f(u) = u *(e^( u)) at (Subscript[u, 0]) = -1, the inverse function g is Lambert W and y(u) = W(u *(e^( u))) on intervals with right hand enpoint -1. This results in the following limit, which we will use twice in what follows.

Corollary.  \(TraditionalForm\`lim\_\(u -> \(-\(1\^-\)\)\)\(W(u*e\^u) + 1\)\/\(u + 1\) = \(-1\)\).

Returning to our projectile with resistance, we now seek a more analytical check of our work by investigating the limiting value for our range formula as the resist ance k goes to 0. Again setting u = -1 - ((k v)/g) sin(θ) and noting that (1/k) = - ((v sin(θ))/(g(u + 1))), we have from our earlier range formula 

	R(θ)=\(\(1\/k\) 
                v\ \(cos(θ)\) \((1 - \(W(u*e\^\(\(\ \)\(u\)\))\)\/u)\)\)=\(\(-\(\(\(v\^2\) \(sin(θ)\) \(cos(θ)\)\)\/\(g(1 + u)\)\)\) \((1 - \(W(u*e\^\(\(\ \)\(u\)\))\)\/u)\)\)=((v^2)/g)sin(θ)cos(θ)((W(u*(e^( u)))-u)/(u*(u + 1))).


Since u approaches -1 from below as k approaches 0 from above, the limit we seek to evaluate becomes

	\(lim\_\(u -> \(-\(1\^-\)\)\)R(θ)\)=(((v^2))/g)sin(θ)cos(θ)* \(lim\_\(u -> \(-\(1\^-\)\)\)\)((W(u*(e^( u)))-u)/(u(u + 1))) .
	
We now show that the limit on the right hand side is 2, so that our results are fully consistent with the no resistance case discussed earlier. Indeed,
	
\(lim\_\(u -> \(-\(1\^-\)\)\)\)((W(u*(e^( u)))-u)/(u(u + 1)))=\(lim\_\(u -> \(-\(1\^-\)\)\)\)((W(u*(e^( u)))-(-1)-(u + 1))/(u(u + 1)))=\(lim\_\(u -> \(-\(1\^-\)\)\)\)((W(u*(e^( u)))-(-1))/(u(u + 1)))+1=-\(lim\_\(u -> \(-\(1\^-\)\)\)\)((W(u*(e^( u)))-(-1))/(u - (-1)))+1.

Applying the Corollary to Theorem 1, the final limit has a value of -1and an overall limit of 2 as claimed. The result is the no resistance range formula R(θ)=(((v^2))/g)*2sin(θ)cos(θ)=((v^2)/g)sin(2  θ).	
	
It is worth noting that Mathematica seems to offer a much more direct proof as it symbolically evaluates the limit we started with. The first command that follows defines Lambert's W function, built in to Mathematica as the function ProductLog. The second command then evaluates, symbollically, the desired one-sided limit.

W = ProductLog;
Limit[((W[u*Exp[u]] - u)/(u (u + 1))), u -> -1, 
      Direction -> 1]

The authors are not aware of precisely what method Mathematica uses to evaluate this limit and we have not found an "easy" proof (this was our motivaton for Theorem 1). Despite our general enthusiasm for machine computation, we would not be comfortable with publication use of a "Mathematica proof" at this time. It is interesting to speculate, however, on how mathematical standards for proof might change if and when "experimental mathematics" and symbolic computation become more widely appreciated and accepted.

The Inverse Range Problem

In the spirit of Groetsch [4], we now consider the inverse problem of what can be proved about angle(s) θ that satisfy R(θ) = r for a given range value r. While finding a closed form inverse relationship does not seem feasible at this point (perhaps it awaits the invention of a few more functions), we can use the form of R(θ) to prove much of what is apparent from our graphs and our physical intuition. The notation and the lemma that follow will streamline our efforts in this regard, which involve extensive use of the chain rule and some surprisingly compact and helpful identities. 

Notation.
 u   =  -1 - (v/g) 
            k  sin (θ)  
 y   =  W(u (e^u))  
 R   =  
        v cos( θ ) (1/k) (1 - (y/u)), where we treat k and θ as the independent variables.
 
Lemma. 
 
(a)  -1 < y < 0 and 0 < (y/u) < 1 for all 0 < θ < (π/2) and hence for all u < -1.

(b)  ((d y)/(d u))=(((u + 1) y)/((y + 1) u))=(((1 + (1/u)))/((1 + (1/y)))) < 0 for all u < -1.

(c)  For u < -1, \(TraditionalForm\`\(\(\(\(d\)\(\ \)\)\^2\) y\)\/\(d\ u\^\(\(\ \)\(2\)\)\) < 0\) and -1<((d y)/(d u)) < 0.

(d) (d/(d u)) (1 - 
              (y/u))  = - (1/(u + 1)) ((d y)/(d u)) (1 - 
              (y/u))

(e)  (∂_k) 
          R   = -(1/(k^2)) v cos( θ) (1 - 
              (y/u)) (1 + ((d y)/(d u)))

(f)    \(TraditionalForm\`∂\_θ 
          R\ \  = \(\(-\(v\/k\)\) 
            sin\ \((\ θ)\)\ \((1 - 
                y\/u)\)\ \((1 + \ \(d\ y\)\/\(d\ u\)\ \(\(cos\^2\)(\ θ)\)\/\(sin\^2\ \((\ θ)\)\))\) = \(-\(v\/k\)\) 
            sin\ \((\ θ)\)\ \(\(cot\^2\)(\ θ)\)\ \((1 - 
                y\/u)\)\ \(\((\(tan\^2\)(\ θ) + \ \(d\ y\)\/\(d\ u\))\)\(.\)\)\)\)


Proofs.

(a)  Since -(1/e) < u (e^u)<0 for u < -1, we get -1 < W(u (e^u)) = y < 0 and consequently 0 < (y/u) < 1.

(b)  Since W is the inverse to T(u) = u (e^u) on [-1,  ∞), for any u < -1we get  T(y) = T(W(u (e^u))), which says y (e^y) = u (e^u). Implicit differentiation then gives ((d y)/(d u))=(((u + 1) (e^u))/((y + 1) (e^y)))=(((u + 1) y)/((y + 1) u)). Using u < -1 and (a), we see that ((d y)/(d u)) has 3 negative factors and a positive factor and must be negative.

(c)  Note that by the Corollary to Theorem 1, \(TraditionalForm\`lim\+\(x -> \(-\(1\^-\)\)\)\(y + 1\)\/\(u + 1\)\  = \ \(-1\)\).  Using (b),  \(TraditionalForm\`lim\+\(u -> \(-\(1\^-\)\)\)\(d\ y\)\/\(d\ u\) = \(lim\+\(u -> \(-\(1\^-\)\)\)\((\(u + 1\)\/\(y + 1\))\) \(lim\+\(x -> \(-\(1\^-\)\)\)\((y\/u)\)\)\  = \ \(-1\)\)\). Now consider the second derivative, which we also obtain by implicit differentiation ( of \(TraditionalForm\`\(d\ y\)\/\(d\ u\) = \((1 + 1\/x)\)\/\((1 + 1\/y)\)\)) and substitution back in for ((d y)/(d u)), \(\(\(\(\(d\)\(\ \)\)\^2\) y\)\/\(d\ u\^\(\(\ \)\(2\)\)\) = \(\(\(-\(1\/u\^2\)\) \((1 + 1\/y)\) + \(1\/y\^2\) \((1 + 1\/u)\) y'\)\/\((1 + 1\/y)\)\^2 = \(\(\(-\(1\/u\^2\)\) \((1 + 1\/y)\)\^2 + \(1\/y\^2\) \((1 + 1\/u)\)\^2\)\/\((1 + 1\/y)\)\^3 = \(\((u + 1)\)\^2 - \((y + 1)\)\^2\)\/\(\(u\^2\) \(\(y\^2\)(1 + 1\/y)\)\^3\)\)\)\). Note that the denominator is negative for u < -1.  To prove that the numerator is always positive, let G(u) = (((u + 1))^2) - (((y + 1))^2).  Then G' (u) = 2 (u + 1) - 
            2 (y + 1) y' = 2 (u + 1) - 
              2 (y + 1) (u + 1) (y/((y + 1) u)) = 
            2 (u + 1) (1 - (y/u)). Then G' (u) is easily seen to be negative for u < -1.   Since G(-1) = 0, this implies G(u) > 0 for u < -1.  Thus the numerator of \(TraditionalForm\`\(\(\(\(d\)\(\ \)\)\^2\) y\)\/\(d\ u\^\(\(\ \)\(2\)\)\)\) is always positive.  Thus \(TraditionalForm\`\(\(\(\(d\)\(\ \)\)\^2\) y\)\/\(d\ u\^\(\(\ \)\(2\)\)\)\) is always negative, for u < -1.  Finally, since  \(TraditionalForm\`lim\+\(u -> \(-\(1\^-\)\)\)\(d\ y\)\/\(d\ u\)\  = \ \(-1\)\)and \(TraditionalForm\`\(\(\(\(d\)\(\ \)\)\^2\) y\)\/\(d\ u\^\(\(\ \)\(2\)\)\)\)is negative for u < -1, we have ((d y)/(d u)) > -1 for u < -1.

(d)  By product rule, we have \(TraditionalForm\`d\/\(d\ u\)\ \((1 - 
              y\/u)\)\  = \ \(-\ \(1\/u\)\)\ \(d\ y\)\/\(d\ u\) + \(\(y\/u\^2\)\(.\)\)\) From (b), \(TraditionalForm\`\(d\ y\)\/\(d\ u\) = \(\((u + 1)\) y\)\/\(\((y + 1)\) u\)\), so \(TraditionalForm\`y\/u\  = \((y + 1)\)\/\((u + 1)\)\ \(d\ y\)\/\(d\ u\)\). Thus we have \(TraditionalForm\`d\/\(d\ u\)\ \((1 - 
              y\/u)\)\  = \(\(-\ \(1\/u\)\)\ \(d\ y\)\/\(d\ u\) + \(1\/u\) \((y + 1)\)\/\((u + 1)\)\ \(d\ y\)\/\(d\ u\) = \(\(-\ \(1\/u\)\) \(\(d\ y\)\/\(d\ u\)\) \((1 - \(y + 1\)\/\(u + 1\))\) = \(\(-\ \(1\/u\)\) \(\(d\ y\)\/\(d\ u\)\) \((\(u - y\)\/\(u + 1\))\) = \(-\ \(1\/\(u + 
                        1\)\)\) \(\(d\ y\)\/\(d\ u\)\) \(\((1 - 
                      y\/u)\)\(.\)\)\)\)\)\) 

(e)  By product rule and chain rule, we have (∂_k) ( 
            v cos θ (1/k) (1 - (y/u)))  =  
        v cos θ [- (1/(k^2)) (1 - (y/u))  +  
                (1/k) ((d )/(d u)) (1 - 
                      (y/u)) (∂_k) u]. Using (d) and the fact that (∂_k) 
          u  = -(v/g) sin θ = ((u + 
                1)/k) , we have \(TraditionalForm\`∂\_k\ \((\ 
            v\ cos\ θ\ \(1\/k\) \((1 - 
                  y\/u)\))\)\  = \ \(v\ cos\ θ\ [\(-\ \(1\/k\^2\)\) \((1 - y\/u)\)\  - \ \(1\/k\) 
                  1\/\(u + 1\)\ \(d\ y\)\/\(d\ u\)\ \((1 - 
                      y\/u)\)\ \(u + 1\)\/k] = 
          v\ cos\ \(θ\ [\(-\ \(1\/k\^2\)\) \((1 - y\/u)\)\ ]\)[
              1 + \ \(d\ y\)\/\(d\ u\)]\)\) as claimed.

(f)  By product rule and chain rule, we have (∂_θ) ( 
            v cos θ (1/k) (1 - 
                  (y/u)))  =  (v/k) [- 
                  sin θ (1 - (y/u)) +  
              cos θ ((d )/(d u)) (1 - 
                    (y/u)) (∂_θ) u]. Using (d) and the fact that (∂_θ) 
          u  = -(v/g) k cos θ  =  (u + 
                1) ((cos θ)/(sin θ)) , we have \(TraditionalForm\`∂\_θ\ \((\ 
            v\ cos\ θ\ \(1\/k\) \((1 - 
                  y\/u)\))\)\  = \ \(\(v\/k\)\ [\(-\ 
                  sin\)\ θ\ \((1 - y\/u)\) + \ 
              cos\ θ\ \(-1\)\/\(u + 1\)\ \(d\ y\)\/\(d\ u\)\ \((1 - 
                    y\/u)\)\ \((u + 
                    1)\)\ \(cos\ θ\)\/\(sin\ θ\)] = \(-\(v\/k\)\) sin\ θ\ \((1 - y\/u)\)\ [
              1 + \ \(d\ y\)\/\(d\ u\)\ \(\(cos\^2\) θ\)\/\(sin\^2\ θ\)]\)\) and the second form follows directly.

Theorem 2. For fixed angle θ the range values decrease as the resistance k increases. Also, for any resistance k > 0 there is a maximum range value, M, that occurs for a unique angle (Subscript[θ, max]) < (π/4). Finally, for all nonnegative ranges r < M, there are exactly two angles θ ∈ [0, (π/2)] that satisfy R(θ) = r, one on either side of (Subscript[θ, max]). 

Proof.  By parts (e) and (a) of the lemma, the partial of R with respect to k is negative, so the range decreases with increasing k. Looking at the ((tan^2)( θ) +  ((d y)/(d u))) factor of (∂_θ) R given in part (f) of the lemma, we note that  ((d y)/(d u)) is an increasing function of θ by (b) of the lemma ((((d^2) y)/(d (u^2))) < 0) and the fact that u is a decreasing function of θ. Since (tan^2)is also increasing in θ, \(TraditionalForm\`\((\(tan\^2\) θ + \ \(d\ y\)\/\(d\ u\))\)\) is an increasing function of θ. Since \(TraditionalForm\`\((\(tan\^2\) θ + \ \(d\ y\)\/\(d\ u\))\)\) goes to -1 as \(TraditionalForm\`θ -> \(0\^+\)\)and goes to ∞ as \(TraditionalForm\`θ -> \(\((π\/2)\)\^-\)\), we conclude that ((tan^2) θ +  ((d y)/(d u))) = 
        0 at a unique (Subscript[θ, 0]) strictly between 0 and (π/2). Noting that for 0 < θ < (π/2), we have (∂_θ) R = 0 if and only if ((tan^2) θ +  ((d y)/(d u))) = 
        0 and we use -1 < ((d y)/(d u)) < 0 (part (c) of the lemma) to conclude that the unique (Subscript[θ, max]) where  \(TraditionalForm\`\((\(tan\^2\) θ\_max + \ \(d\ y\)\/\(d\ u\))\) = 0\) must satisfy (Subscript[θ, max]) < (π/4). The fact that this unique (Subscript[θ, max]) maximizes R follows since (∂_θ) R is positive to the left of (Subscript[θ, max]) and negative to its right. That there are exactly two elevation angles that achieve any submaximal range also follows from what was just said about the signs of (∂_θ) R and the fact that R is zero when θ = 0 and when θ = (π/2).

A Closed Form for the Maximum Range Angle

The efforts of the previous section provide a surprising additional dividend. Using some interesting algebra and our old friend Lambert W, we are able to find a formula for the angle that maximizes the range for given values of the resistance k, the initial velocity v, and the acceleration of gravity g. In fact, the key parameter turns out to be α = ((k v)/g), a notation we will use for the remainder of this section. A similar result using a "log-like" function has recently been offered by Groetsch [6]. With some further analysis are able to show that the maximum range function is continuous and decreasing in α. 

Theorem 3.  The angle (Subscript[θ, max]) that gives the maximum range for a given  v, g, k depends only on α = ((k v)/g) and is given by the formula  (Subscript[θ, max])=ArcSin[((αW((((α^2) - 1)/e)))/((α^2) - 1 - 
                  W((((α^2) - 1)/e))))] for α ≠ 1, and (Subscript[θ, max])=ArcSin[(1/(e - 1))] for α = 1 . 

Proof.  For simplicity, just write θ = (Subscript[θ, max]). From Theorem 2, we have \(TraditionalForm\`\(\(sin\^2\) θ\)\/\(\(cos\^2\) θ\) + \ \(d\ y\)\/\(d\ u\) = 0\). Using \(TraditionalForm\`\(d\ y\)\/\(d\ u\) = \(\((u + 1)\) y\)\/\(u(y + 1)\)\) and sin θ = -((u + 1)/α), this can be rewritten as \(TraditionalForm\`\(\(1\/α\^2\) \((u + 1)\)\^2\)\/\(1 - \(1\/α\^2\) \((u + 1)\)\^2\) = \(-\(\(\((u + 1)\) y\)\/\(u(
                y + 1)\)\)\)\).  Then \(TraditionalForm\`\(\(\ \)\(\(u(u + 1)\)\/\(α\^2 - \((u + 1)\)\^2\) = \(-\(y\/\((y + 1)\)\)\)\)\)\).  Taking reciprocals, \(TraditionalForm\`\(\(\ \)\(\(α\^2 - \((u + 1)\)\^2\)\/\(u(u + 1)\) = \(-\(\(y + 1\)\/y\)\)\)\)\).  Adding 1 to both sides, \(TraditionalForm\`\(α\^2 - \((u + 1)\)\^2 + u(u + 1)\)\/\(u(u + 1)\) = \(-\(1\/y\)\)\).  So \(TraditionalForm\`\(α\^2 - \((u + 1)\)\)\/\(u(u + 1)\) = \(-\(1\/y\)\)\) and y = -((u(u + 1))/((α^2) - 1 - 
                u)).  Now, using y (e^y) = u (e^u),  we get -((u(u + 1))/((α^2) - 1 - u)) 
          (e^(-((u(u + 1))/((α^2) - 1 - u)))) = u (e^u).  Then -(((u + 1))/((α^2) - 1 - u)) = 
        (e^(u + ((u(u + 1))/((α^2) - 1 - u)))), and so -(((u + 1))/((α^2) - 1 - u)) = 
        (e^((((α^2) u)/((α^2) - 1 - u)))).  

Case α = 1:  We have ((u + 1)/u) = (e^(-1)).  Solving for u leads to  u=((-e)/(e - 1)).  Then sin (Subscript[θ, max])=-((u + 1)/α)=(1/(e - 1)), as claimed.

Case α ≠ 1:  To solve this for u,  make the substitution  t = (1/((α^2) - 1 - u)) (and so u = (α^2) - 1 - (1/t)) to get -((α^2) - (1/t)) t=(e^( (α^2)((α^2) - 1 - (1/t)) 
              t)).  So (α^2) t-1+(e^( (α^2)((α^2) - 1) 
                  t - (α^2)))=0.  This equation is in the form given in result (3) of the section where the Lambert W function was introduced with  a = (α^2), b = -1, c=(e^( -(α^2))), d = (α^2)((α^2) - 1).  Result (3) therefore gives t=\(1\/α\^2 - \(1\/\(\(α\^2\)(α\^2 - 
                        1)\)\) \(W(\(\(e\^\(-α\^2\)\) \(\(α\^2\)(α\^2 - 1)\) e\^\(\(\(α\^2\)(α\^2 - 1)\)\/α\^2\)\)\/α\^2)\)\)=(1/(α^2))-(1/((α^2)((α^2) - 1)))W((((α^2) - 1)/e)).  Using u = (α^2) - 1 - (1/t) again, we get u=(α^2)-1-(1/((1/(α^2))-(1/((α^2)((α^2) - 1)))W((((α^2) - 1)/e))))=(α^2)-1-(((α^2)((α^2) - 1))/((α^2)-1-W((((α^2) - 1)/e)))).  Then sin θ=-((u + 1)/α)=-(1/α)((α^2)-(((α^2)((α^2) - 1))/((α^2)-1-W((((α^2) - 1)/e))))), which yields the formula sin( θ)=((α W((((α^2) - 1)/e)))/((α^2)-1-W((((α^2) - 1)/e)))), as desired.


In Figure 6 we show a plot of the maximum angle function as a function of α. In the spirit of this paper, we will back up this visual evidence of technology with some careful analysis.

Figure 6.  Maximum elevation angle as a function of α = ((k v)/g)

We close this section by establishing two properties of the (Subscript[θ, max]) function. The first is the physically nonsurprising fact that this function is continuous for all α ≥ 0. Inspecting the formula for (Subscript[θ, max]) from Theorem 3, it suffices to show that  (Subscript[θ, max]) is continuous at α = 0 and α = 1.  First we need to show that \(TraditionalForm\`lim\+\(α -> \(0\^+\)\)θ\_max = π\/4\). The key is that if we write W=W((((α^2) - 1)/e)), then W (e^W)=(((α^2) - 1)/e) (using result (1) for Lambert W).  Now expand the left side as a power series about W = -1.  (This power series converges for all real numbers W.)  So -(1/e)+(1/(2e))(((W + 1))^2)+(2/(6e))(((W + 1))^3)+...=(((α^2) - 1)/e). Thus \(TraditionalForm\`\(1\/2\) \((W + 1)\)\^2 + \(2\/6\) \((W + 1)\)\^3 +  ..  = α\^2\).  If we denote  H(W) = (1/2) + (2/6) (W + 1) +  ... , then  H(W) also converges for all real numbers, and in particular is continuous at W = -1. Note (α^2) = (((W + 1))^2) H(W).  From this, we get \(TraditionalForm\`lim\+\(α -> \(0\^+\)\)\(W + 1\)\/α = \(\@\(lim\+\(α -> \(0\^+\)\)\((W + 1)\)\^2\/α\^2\) = \(\@\(lim\+\(W -> \(-\(1\^+\)\)\)\((W + 1)\)\^2\/\(\(\((W + 1)\)\^2\) \(\(H(W)\)\(.\)\)\)\) = \(\@\(lim\+\(W -> \(-\(1\^+\)\)\)1\/\(\(H(W)\)\(.\)\)\) = \@2\)\)\)\). Then we have \(lim\+\(α -> \(0\^+\)\)sin(\ θ\_max)\)=\(lim\+\(α -> \(0\^+\)\)\)((αW((((α^2) - 1)/e)))/((α^2)-1-W((((α^2) - 1)/e))))=\(lim\+\(α -> \(0\^+\)\)\)((W((((α^2) - 1)/e)))/(α-((1+W((((α^2) - 1)/e)))/α)))=((-1)/(0 - (Sqrt[2]))) = (1/(Sqrt[2])).  This implies \(TraditionalForm\`lim\+\(α -> \(0\^+\)\)θ\_max = π\/4\), which is what we expect as the resistance k -> 0.  We must also show that (lim_(α -> 1)) ((α W((((α^2) - 1)/e)))/((α^2)-1-W((((α^2) - 1)/e))))=(1/(e-1)) .   This can be done by rewriting the limit as (lim_(α -> 1)) (α/((((α^2) - 1)/(W((((α^2) - 1)/e))))-1)) and using L'Hopitals Rule on the part (lim_(α -> 1)) (((α^2) - 1)/(W((((α^2) - 1)/e))))=(lim_(α -> 1)) ((2  α)/(W'((((α^2) - 1)/e))2(α/e)))=(e/(W' (0)))=(e/1)=e. Thus (lim_(α -> 1)) (α/((((α^2) - 1)/(W((((α^2) - 1)/e))))-1))=(1/(e-1)), as desired. We note in passing that Mathematica had no trouble (other than being "slow") in computing \(TraditionalForm\`lim\+\(α -> \(0\^+\)\)sin(\ θ\_max) = 1\/\@2\), but was unable to evaluate \(TraditionalForm\`lim\+\(α -> 1\)\(\(sin(\ θ\_max)\)\(.\)\)\)


The second property of (Subscript[θ, max]) confirms our physical intuition that increasing either the resistance or the initial velocity will decrease the maximum range angle. More precisely (and pleasingly), (Subscript[θ, max]) is a decreasing function of α. To show this, it is convenient to look at the reciprocal \(TraditionalForm\`1\/\(sin(\ θ\_max)\) = \(α\^2 - 1\)\/\(α\ W\) - 1\/α\). The derivative of this (again using θ for (Subscript[θ, max])) is 
(d/(d α)) ((1/(sin( θ))))=(1/((α^2) (W^2)))[2  α α W-((α^2) - 1)(W+α(W/((((α^2) - 1)/e)(1 + W)))2(α/e))]+(1/(α^2))=(1/((α^2) (W^2)(1 + W)))[
                2  (α^2) (W^2) + W + 1]. Note that we used the fact ((d W)/(d α))=((d )/(d α))W((((α^2) - 1)/e))=(W/(((((α^2) - 1)/e))(W + 1)))2(α/e) by the Chain Rule. The above expression for \(TraditionalForm\`\(d\/\(d\ α\)\) \((1\/\(sin(\ θ)\))\)\) is clearly positive for α > 0 and α ≠ 1 (because W + 1 > 0) . Since we already know (1/(sin( θ))) is a continuous function on α ≥ 0, then (1/(sin( θ))) is an increasing function of α.  Since (1/(sin( θ)))  is also positive function, then sin(θ) is a decreasing function of α. Since θ is between 0 and (π/2),  (Subscript[θ, max]) is therefore a decreasing function of α.

Concluding Comments

Earlier we mentioned the phrase Experimental Mathematics, which is relevant to this paper in several ways. It was only in working with Mathematica that we realized that the Lambert W function offered a way of getting a closed form solution for the range function in the presence of linear resistance. Also, symbolic and graphical outputs suggested some of our formal mathematical efforts and swiftly confirmed others. In our belief that this may become an increasingly useful approach for doing mathematical research, we quote part of the statement of philosophy of the Journal of Experimental Mathematics, a journal started in 1992.

Experiment has always been, and increasingly is, an important method of mathematical discovery. ....  Yet this tends to be concealed by the tradition of presenting only elegant, well-rounded and rigorous results. 

While we value the theorem-proof method of exposition, and while we do not depart from the established view that a result can only become part of mathematical knowledge once it is supported by a logical proof, we consider it anomalous that an important component of the process of mathematical creation is hidden from public discussion. It is to our loss that most of us in the mathematical community are almost always unaware of how new results have been discovered. ....

Experimental Mathematics was founded in the belief that theory and experiment feed on each other, and that the mathematical community stands to benefit from a more complete exposure to the experimental process. The early sharing of insights increases the possibility that they will lead to theorems: an interesting conjecture is often formulated by a researcher who lacks the techniques to formalize a proof, while those who have the techniques at their fingertips have been looking elsewhere. Even when the person who had the initial insight goes on to find a proof, a discussion of the heuristic process can be of help, or at least of interest, to other researchers. There is value not only in the discovery itself, but also in the road that leads to it.

One might question the need for detailed and sometimes complex proofs of results which which are either physically "obvious" or obtainable by a simple sequence of commands in a computer algebra system. We would respond by not only reaffirming the importance of careful proof, but also by emphasizing the thrill of the chase and the satisfaction when results work out so cleanly and lead to more general results. In particular, Theorem 1 evolved from the specific need to check results of our "experimental" computer investigation. And we were both surprised and pleased at the concise final form of Theorem 3 and the crisp identities that came up in the process of its proof.

The success of the Lambert W function in obtaining closed form solutions suggests several directions for further projectile research. We have already mentioned, without much optimism, the inverse problem of finding what θ value(s) give rise to a given range R. One might also consider the more general overall problem of firing a projectile at a target located on an inclined plane (Groetsch, [5]) or the sad but likely fact that resistance contributes a term that is probably not linear in terms of the velocity. Regardless of its possible help in attacking these and other projectile problems, we have become believers in Lambert W and are confident that its role in mathematics will continue to expand.

References

1.  ---, Time for a New Elementary Function?, FOCUS; Newsletter Math. Assoc. Amer. 20 (2000) 2. 
2.  Jonathan M. Borwein and Robert M. Corless, Emerging Tools for Experimental Mathematics, Amer. Math Monthly 106 (1999) 889-909.
3.  R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, On the Lambert W Function, Adv. Comput. Math. 5 (1996) 329-359.
4.  C. W. Groetsch, Tartaglia's Inverse Problem in a Resistive Medium, Amer. Math. Monthly 103 (1996) 546-551.
5.  C. W. Groetsch, Halley's Gunnery Rule, The College Math. Journal 28 (1997) 49-50.
6.  C. W. Groetsch, Tartaglia's Bet, forthcoming in Cubo Matematica Educacional.