Goal

The goal of this lab is to use the Euler method to model system that are changing with time. If we can describe how the system changes with time, then we can determine how the system changes by stepping it forward in time using tiny intervals of time. This is the basic idea behind the Euler method. We apply this method to two-dimensional motion.

Overview

In order to visualize the motion of a cannonball projectile with effects of air resistance, we utilize the Euler method to obtain position and velocity in coordinate directions. The steps necessary to do this has been compartmentalized as shown below:

• Initial Values
• Acceleration in $x$ and $y$
• Calculating position and velocity
• Plotting the results

Setting initial values

To begin our problem, we set initial values for the variables we need in our calculation. This includes initial position and velocity as well as various constants.

# enter the mass, radius of cannon ball
m = 20 # what is a good value for the mass of a cannon ball
r = 0.3 #radius of cannon ball

# enter acceleration
g = -9.81

# enter initial conditions
theta0_deg = 30 # degrees
v0 = 700. #m/s for a cannon ball
# convert angle to radians - numpy functions need angle in radians!
theta0 = np.radians(theta0_deg)
y0 = 100
x0 = 0

# calculate the initial x and y velocities
v0x = v0*np.cos(theta0)
v0y = v0*np.sin(theta0)

# enter time step dt
dt = 0.01

For air resistance, $B_2$ is given by:

$$B_2 = \frac{1}{2}C\rho A$$

where $C$ is a constant, $\rho$ is the density of air, and $A$ is cross-sectional area.

 # calculate B2 here and save it as myB2
C = 0.5
rho = 1.225
A = np.pi*r**2
myB2 = 0.5*C*rho*A

Calculating acceleration

In order to calculate acceleration in $x$ and $y$, we define two functions to do so. $a_{drag}$ is given by:

$$a_{drag,x,y} = -\frac{B_2}{m}v v_{x,y}$$

$a_{drag,y}$ also includes gravitational acceleration: $g$.

# write your function to calculate the horizontal acceleration
def accel_x(v,vx,m,myB2):
  '''calculates acceleration in the x direction'''
  ax = -myB2/m*v*vx
  return ax

# write your function to calculate the vertical acceleration
def accel_y(v,vy,m,myB2):
  '''calculates acceleration in the y direction'''
  ay = g-myB2/m*v*vy
  return ay

Calculating position and velocity

To calculate the position and velocity of the cannonball for each time step, we define a function to calculate these values.

# write your calculate function to compute the Euler solution for position and velocity
# NOTE: this function will initialize the x/y position and velocities as empty lists
def calculate(x0,y0,v0x,v0y,v0,dt):
  '''calculates x,y,vx,vy using intial conditions'''
  x = [x0]
  y = [y0]
  vx = [v0x]
  vy = [v0y]
  i = 1
  v = [v0]
  while (y[i-1]>=0):
    x.append(x[i-1] + vx[i-1]*dt)
    y.append(y[i-1] + vy[i-1]*dt)
    vx.append(vx[i-1] + accel_x(v[i-1],vx[i-1],m,myB2)*dt)
    vy.append(vy[i-1] + accel_y(v[i-1],vy[i-1],m,myB2)*dt)
    v.append(np.sqrt(vx[i-1]**2+vy[i-1]**2))
    i=i+1
  return x,y,vx,vy

We used the while statement to limit the range of our function for values of $y$ at and above 0 (cannonball's position above ground).

Plotting

After we run our calculate function, we can plot the resulting values to view the trajectory. These were the resulting plots from the initial conditions set above:

We can also use our function to view the effects of differing launch angles:

#calculates thetas
theta0_deg = np.arange(15,90,15)
theta0 = np.radians(theta0_deg)
#calculates initial velocities
v0x = v0*np.cos(theta0)
v0y = v0*np.sin(theta0)

x1,y1,vx1,vy1 = calculate(x0,y0,v0x[0],v0y[0],v0,dt)
x2,y2,vx2,vy2 = calculate(x0,y0,v0x[1],v0y[1],v0,dt)
x3,y3,vx3,vy3 = calculate(x0,y0,v0x[2],v0y[2],v0,dt)
x4,y4,vx4,vy4 = calculate(x0,y0,v0x[3],v0y[3],v0,dt)
x5,y5,vx5,vy5 = calculate(x0,y0,v0x[4],v0y[4],v0,dt)

plt.figure()
plt.plot(x1,y1,label="15 deg")
plt.plot(x2,y2,label="30 deg")
plt.plot(x3,y3,label="45 deg")
plt.plot(x4,y4,label="60 deg")
plt.plot(x5,y5,label="75 deg")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()

Resultant plot:

Without air resistance (set myB2=0 and use above code):

Note that with air resistance, lower launch angle results in greater range. This is unlike without air resistance, where 45 degrees obtains the furthest range.