Assignment 4
General Feedback
Some students solved the “no air resistance” problem analytically. It is more practical to use the existing code in the relevant limit (kappa = 0 or rho_air = 0).
Some students did not read the plotting instructions carefully and plotted the trajectories for a second time but this time without drag.
Some students were still missing labels, title or a legend
Model solution
Automatically marked tasks
def get_area(r):
'''
This function returns the cross section area of a sphere of radius r. The returned
value is in the squared units of the unit of the radius argument.
'''
### BEGIN SOLUTION ###
area = numpy.pi * r**2 # Cross sectional area of sphere
return area
### END SOLUTION ###
def get_mass(r):
'''
This function returns the mass of an iron sphere of radius r. The radius
should be given in meter and the return value is in kg.
'''
### BEGIN SOLUTION ###
mass = rho_iron * (4/3) * numpy.pi * r**3 # Mass of sphere
return mass
### END SOLUTION ###
area_cb = get_area(r_cb)
mass_cb = get_mass(r_cb)
def f(r, t):
'''Implements differential equation for cannonball from state vector r and time t'''
# Unpack array of the state
x, y, vx, vy = r
# these variables should updated in your code to be the derivatives of
# the x, y positions and the derivative of the x, y velocities.
dx_t, dy_dt, dvx_dt, dvy_dt = 0, 0, 0, 0
### BEGIN SOLUTION ###
v = numpy.sqrt(vx**2 + vy**2)
k = 0.5 * rho_air * kappa * area_cb
# x-components
dx_dt = vx # DE for x position
dvx_dt = -(k/mass_cb)*v*vx # DE for x velocity
# y-components
dy_dt = vy # DE for y position
dvy_dt = -(k/mass_cb)*v*vy - g
### END SOLUTION ###
return numpy.array([dx_dt, dy_dt, dvx_dt, dvy_dt])
def solve_euler(state_initial, t1, n_steps):
'''Solves ODE using Euler's method from state_initial to end time t1 using n_panels panels'''
# Define array for trajectory history
history = numpy.empty((n_steps+1,4)) # Number of columns is equal to number of variables to solve for (4)
history[0] = state_initial
# you should now populate the history array
### BEGIN SOLUTION ###
# Define necessary local variables
t = 0
dt = t1 / n_steps # Panel size
for i in range(n_steps):
t = i*dt
history[i+1] = history[i] + f(history[i], t) * dt
### END SOLUTION ###
return history
def find_zero_linear(x1, x2, y1, y2):
if y1*y2 > 0:
print("I expect y1 and y2 to have opposite signs!")
### BEGIN SOLUTION ###
return (x1*y2 - x2*y1 )/(y2 - y1)
### END SOLUTION ###
Plotting tasks
These are examples for the plotting tasks.
n_steps = 1000
max_time = 300
v0s = numpy.linspace(50, 1000, 20)
ranges = []
ranges_noresistance = []
theta = numpy.deg2rad(60)
v0 = 125
n_steps = 1000
thetas = range(5, 90, 5)
t1=100
# Loop creating trajectories
for theta in thetas:
theta = numpy.deg2rad(theta) # numpy trig functions work in radians
vx, vy = v0*numpy.cos(theta), v0*numpy.sin(theta)
initial_conditions = [0, 0, vx, vy]
values_euler = solve_euler(initial_conditions, t1, n_steps)
# Plot trajectories
xs_euler, ys_euler = values_euler[:,0], values_euler[:,1]
plt.plot(xs_euler, ys_euler, color='blue', linestyle='--')
plt.xlim(0,1500)
plt.ylim(0,800)
plt.title('Trajectories at Different Launch Angles')
plt.xlabel('$x$ (m)')
plt.ylabel('$y$ (m)')
n_steps = 1000
max_time = 300
v0s = numpy.linspace(50, 1000, 20)
ranges = []
ranges_noresistance = []
theta = numpy.deg2rad(60)
for v0 in v0s:
vx, vy = v0*numpy.cos(theta), v0*numpy.sin(theta)
initial_conditions = [0, 0, vx, vy]
rho_air = 0
history = solve_euler(initial_conditions, max_time, n_steps)
rng = find_range(history)
ranges_noresistance.append(rng)
rho_air = 1.23
history = solve_euler(initial_conditions, max_time, n_steps)
rng = find_range(history)
ranges.append(rng)
plt.plot(v0s,ranges_noresistance,label="No resistance")
plt.plot(v0s,ranges, label="with resistance")
plt.title('Range as a function of the initial velocity')
plt.xlabel('$v_0$ (m/s)')
plt.ylabel('range (m)')
plt.legend()
Using scipy
As with last week’s assignment we can use scipy to solve the problem, it’s a bit more tricky as we need the solution at intermediate time steps.
import scipy.integrate as integrate
def f_swap(t,n) :
return f(n,t)
n_steps = 1000
max_time = 300
v0s = numpy.linspace(50, 1000, 20)
ranges = []
ranges_noresistance = []
theta = numpy.deg2rad(60)
times=numpy.linspace(0,max_time,n_steps+1)
for v0 in v0s:
vx, vy = v0*numpy.cos(theta), v0*numpy.sin(theta)
initial_conditions = numpy.array([0, 0, vx, vy])
rho_air = 0
output = integrate.solve_ivp(f_swap,(0,max_time),initial_conditions,t_eval=times,rtol=1e-9)
history = numpy.transpose(output.y)
rng = find_range(history)
ranges_noresistance.append(rng)
rho_air = 1.23
output = integrate.solve_ivp(f_swap,(0,max_time),initial_conditions,t_eval=times,rtol=1e-9)
history = numpy.transpose(output.y)
rng = find_range(history)
ranges.append(rng)
plt.plot(v0s,ranges_noresistance,label="No resistance")
plt.plot(v0s,ranges, label="with resistance")
plt.title('Range as a function of the initial velocity')
plt.xlabel('$v_0$ (m/s)')
plt.ylabel('range (m)')
plt.legend();
