Projectile Motion Simulation with Aerodynamics using Euler's Method

Question

I'm trying to simulate the motion of a projectile taking into account aerodynamic drag. The code works perfectly if the aerodynamic drag is zero. However, if the drag coefficient or velocity is too high, the trajectory starts to bend towards the negative x axis. I used this source to help me create the simulation.

Here's the code:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

cd = 0.1
v = 150
ang = 40 # in degrees
t = [0]  # list to keep track of time
vx = [v * np.cos(ang / 180 * np.pi)] 
vy = [v * np.sin(ang / 180 * np.pi)]
x = [0]  # list for x and y position
y = [0]
g = 9.8
M = 1
density = 1.2
area = 1
# Drag force
drag = area * density * cd * 0.5 * v ** 2
# Acceleration x and y
ax = [-(drag * np.cos(ang / 180 * np.pi)) / M]
ay = [-g - (drag * np.sin(ang / 180 * np.pi) / M)]
# time-step
dt = 0.01
f1 = plt.figure()
ax1 = f1.add_subplot(111)
line, = ax1.plot([], [], '-')
def animate(i):
    t.append(t[i] + dt) 
    vx.append(vx[i] + dt * ax[i])  # Update the velocity
    vy.append(vy[i] + dt * ay[i])
    # Update position
    x.append(x[i] + dt * vx[i])
    y.append(y[i] + dt * vy[i])
    # Calculate updated velocity
    vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2)  # magnitude of velocity
    drag = area * density * 0.5 * cd * vel ** 2
    ax.append(-(drag * np.cos(ang / 180 * np.pi)) / M)
    ay.append(-g - (drag * np.sin(ang / 180 * np.pi) / M))
    line.set_xdata(x[:-1])
    line.set_ydata(y[:-1])
    return line,

ax1.set_ylabel("y (meters)")
ax1.set_xlabel("x (meters)")
ax1.set_xlim(0, 90)
ax1.set_ylim(0, 90)
ax1.autoscale(False)
ani = animation.FuncAnimation(f1, animate, interval=0.1, blit=True)
plt.show()

This is the output:

It looks pretty accurate up until the peak of the trajectory, where eventually it starts bending inward which is obviously incorrect. Why is this happening? Is there an error in the script, or is it simply the result of Euler's method being inaccurate?

Answer 1

You are not updating the angle in your simulation as the velocity vector is changing. Here is a modified version of your 'animate' function and the result.

def animate(i):
    t.append(t[i] + dt) 
    vx.append(vx[i] + dt * ax[i])  # Update the velocity
    vy.append(vy[i] + dt * ay[i])
    vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2)  # magnitude of velocity
    cos = vx[-1]/vel # note: should deal with vel=0 properly
    sin = vy[-1]/vel
    # Update position
    x.append(x[i] + dt * vx[i])
    y.append(y[i] + dt * vy[i])
    # Calculate updated velocity
    drag = area * density * 0.5 * cd * vel ** 2
    ax.append(-(drag * cos) / M)
    ay.append(-g - (drag * sin / M))
    line.set_xdata(x[:-1])
    line.set_ydata(y[:-1])
    return line,

and the result: