Projectile with Drag

Here is another cool solution, that simulates the trajectory numerically.

Firstly, the differential equations that describe the trajectory can be written in terms of Newton's laws:

$\displaystyle \bold{ \vec{F}} = \bold{\vec{F}}_g + \bold{\vec{F}}_d$ , where $\bold{\vec{F}}_g$ and $\bold{\vec{F}}_d$ are the forces of gravity and drag respectively.

The equation can be written as a second-order ODE in vector form:

$\displaystyle m \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} = - \frac{mg}{v_o} \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} - m \begin{bmatrix} 0 \\ g \end{bmatrix}$

Solving through x and y with the forward Euler method, and trialing with different angles, I found that 38 degrees gave the greatest value of $x$ at peak height.

Here's the code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

import math

time = 0
deltaT = 10**-4 #timestep

x = 0
y = 0

#play around with the value of theta. You'll find that around 38 degrees will give maximum x value.

theta = 38

#initial value of velocity

xDot = 100*math.cos(theta*(math.pi/180))
yDot = 100*math.sin(theta*(math.pi/180))


while yDot >= 0: #when vertical velocity is zero the height is maximum


  #differential equations (values of acceleration in x and y)
  xDotDot = -(10/100)*xDot
  yDotDot = -(10/100)*yDot - 10


  #Euler method
  xDot += xDotDot*deltaT
  yDot += yDotDot*deltaT

  x += xDot*deltaT
  y += yDot*deltaT

  time += deltaT


print("x at maximum height for the value theta = " +str(theta) + " is: " + str(x))

We can write the equations of motion along the horizontal and the vertical directions as $\dfrac{dv_x}{dt}=-αv_x, \dfrac{dv_y}{dt}=-α(v_y+\dfrac{g}{α})$ . Solving these two we get $x=\dfrac{v_0\cos \theta_0}{α}(1-e^{-αt}), v_y=-\dfrac{g}{α}+\left (v_0\sin \theta_0+\dfrac{g}{α}\right )e^{-αt}$ . At the highest point of the path, $v_y=0\implies e^{-αt}=\dfrac{g}{g+αv_0\sin \theta_0}$ . Hence, at this point, $x=\dfrac{v_0\cos \theta_0}{α}\left (1-\dfrac{g}{g+αv_0\sin \theta_0}\right )=\dfrac{v_0^2}{g}\times \dfrac{\sin \theta_0\cos \theta_0}{1+\sin \theta_0}$ (since it is given that $αv_0=g$ ). When $x$ is maximum, $\dfrac{dx}{d\theta_0}=0\implies \sin \theta_0\approx 0.61803\implies \theta_0\approx \boxed {38.1724\degree}$ .