In my previous post, I discussed how to simulate vertical motion of the free fall of a projectile without considering air resistance. In this section, I am going to simulate another kind of projectile motion cosidering effect of air resistance. Consider a cannonball fired towards a target with a certain velocity, say $100 m/s$. Suppose, the target is at a distance 100 kilometer from the cannon. What should be the angle of fire so that cannonball hits the target accurately? Consider, the air resistance to be approximately 0.01.

The projectile motion is influenced by force of gravity and drag force due to air resistance. The equation of motion -

For motion in horizontal direction

$$\frac{d^2x}{dt^2}=-kv_x$$

For motion in vertical direction $$\frac{d^2y}{dt^2}=-g -kv_y$$

For numerically solving the differential equations, let us split them into two first order differential equations.

$$\frac{dx}{dt}=v_x $$ $$\frac{dv_x}{dt}=-kv_x$$ $$\frac{dy}{dt}=v_y $$ $$\frac{dv_y}{dt}=-g-kv_y$$

The initial conditions are $$ x(0) = 0\\ y(0) = 0\\ v_x(0)=v\cos\theta\\ v_y(0)=v\sin\theta $$

$\theta$ is the angle of projection.

So, the plan of simulation is as follows -

Solving the differential equations by any standard method. I am using odeint function from SciPy package to solve the differential equations. For a particular angle of projection ($\theta$), the solver function will give horizontal distance ($x$), vertical distance ($y$), horizontal velocity ($v_x$) and vertical velocity ($v_y$) at different time step ($t$).
By applying linear interpolation, we can exactly estimate the point where the projectile hits the ground.
Then, we can find out the difference between the target and the point of hitting on the ground by the projectile.
The difference existing between the target and the point of hitting on the ground by the projectile for an arbitrary angle of projection can be minimised below a tolerance limit by searching appropriate value of angle. This can be done by "Bisection method" usually employed to locate the root of a polynomial function.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import bisect

#basic parametrs
k = 0.01
g = 9.8
v = 100
target = 500

t_init, t_final, step_size = 0, 20, 0.001
t = np.arange(t_init, t_final, step_size)

z = np.zeros([len(t),4])

def model(z, t, params):
    x, y, vx, vy = z
    dx_dt = vx
    dy_dt = vy
    dvx_dt = -k*vx
    dvy_dt = -g-k*vy
    dz_dt = np.array([dx_dt, dy_dt, dvx_dt, dvy_dt])
    return dz_dt

@np.vectorize
def diff(theta):
    theta = np.radians(theta)
    params = [k, g, theta]
    
    x0, y0, vx0, vy0 = 0, 0, v*np.cos(theta), v*np.sin(theta)
    z0 = [x0, y0, vx0, vy0]

    sol = odeint(model, z0, t, args=(params,))
    x, y, vx, vy = sol[:, 0], sol[:, 1], sol[:, 2], sol[:, 3]

    y = y[y>=0]
    x = x[:len(y)]
    vx = vx[:len(y)]
    vy = vy[:len(y)]

    xground = x[-2] + y[-2]*(x[-1]-x[-2])/(y[-1]-y[-2])
    diff = xground - target
    return diff

def plot():
    fig, ax = plt.subplots(figsize=(8, 5))
    # set the x-spine
    ax.spines['left'].set_position('zero')

    # turn off the right spine/ticks
    ax.spines['right'].set_color('none')
    ax.yaxis.tick_left()

    # set the y-spine
    ax.spines['bottom'].set_position('zero')

    # turn off the top spine/ticks
    ax.spines['top'].set_color('none')
    ax.xaxis.tick_bottom()

plot()
theta = np.arange(10, 90, 0.1)
plt.plot(theta, diff(theta), ':r')
plt.xlabel('$\\theta$', fontsize=14)
plt.ylabel('$\Delta x$', fontsize=14)
plt.xlim(0, 90)
plt.show()

$\Delta x$ which represents the difference between target and projectile hitting point is plotted against different angle of projection. It is seen that there are two angles of projection for which projectile can accurately hit the target. The first one lies in the interval $(10^{\circ}, 20^{\circ})$ and the second one in the interval $(70^{\circ}, 80^{\circ})$. Using bisect function from SciPy, we can find the exact value of angle.

angle1 = bisect(diff, 10, 20)
angle2 = bisect(diff, 70, 80)
print('\n Angle1 = %0.2f'%angle1,'\n Angle2 = %0.2f'%angle2)
diff1 = diff(angle1)
diff2 = diff(angle2)
print('\n Difference for angle1 = %0.3f'%diff1,'\n Difference for angle2 = %0.11f'%diff2)

 Angle1 = 15.26 
 Angle2 = 73.15

 Difference for angle1 = -0.085 
 Difference for angle2 = -0.00000000002

After getting the exact value of angle of projection, using projectile function, we can depict the projectile trajectory as shown in figure.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import bisect

def model(z, t, params):
    x, y, vx, vy = z
    dx_dt = vx
    dy_dt = vy
    dvx_dt = -k*vx
    dvy_dt = -g-k*vy
    dz_dt = np.array([dx_dt, dy_dt, dvx_dt, dvy_dt])
    return dz_dt

def projectile(angle):
    theta = np.radians(angle)
    params = [k, g, theta]

    
    x0, y0, vx0, vy0 = 0, 0, v*np.cos(theta), v*np.sin(theta)
    z0 = [x0, y0, vx0, vy0]

    sol = odeint(model, z0, t, args=(params,))
    x, y, vx, vy = sol[:, 0], sol[:, 1], sol[:, 2], sol[:, 3]

    y = y[y>=0]
    x = x[:len(y)]
    vx = vx[:len(y)]
    vy = vy[:len(y)]
    return x, y

plot()
x1, y1 = projectile(angle1)
x2, y2 = projectile(angle2)
plt.plot(x1, y1, ls='--', color='purple', label='$\\theta$ = %d$^\circ}$'%angle1)
plt.plot(x2, y2, ls='--', color='blue', label='$\\theta$ = %d$^\circ}$'%angle2)
plt.plot(500, 0, 'ro', markersize=10)
plt.plot(0, 0, 'ko', markersize=10)
plt.ylim(0, 500)
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.title('Projectile Motion', fontsize=16)
plt.legend(frameon=False)
plt.annotate('Starting point', xy=(0,0), xytext=(50, 100), arrowprops=dict(arrowstyle='->'), fontsize=14)
plt.annotate('Target', xy=(500,0), xytext=(350, 100), arrowprops=dict(arrowstyle='->'), fontsize=14)
plt.show()