Problem 2. “Pumpkin chucking” is a competition event to see which team can shoot a pumpkin as far as possible, usually with a pneumatic cannon. In this problem we’re going to write a computer program to determine the optimal angle to shoot a spherical projectile to maximize distance, taking into account air resistance (aerodynamic drag).
Without any air resistance, the optimal angle to shoot any projectile is 45 degrees, but once air resistance is added, the path of the object cannot be derived analytically. The trajectory with air resistance is shown below. Notice that it is not a parabola.
The force due to air resistance is a vector that points in the opposite direction of the velocity vector and has length b|⃗v|2, where b is a coefficient that depends on the shape of the object and the density of air. For a sphere-like pumpkin, which weighs 4 kg and has a radius of 0.1 m (9 pounds with an 8" diameter), we will take the coefficient b to be 0.008 kg/m.
(a) Write a (continuous) dynamical system for the 2D position of a pumpkin with mass m shot at time t = 0 from position (0, 0) at an initial speed s and angle θ relative to the flat ground. Your system should be of the form dx1 = f1(x1,x2,x3,x4) dt dx2 = f2(x1,x2,x3,x4) dt dx3 = f3(x1,x2,x3,x4) dt dx4 = f4(x1,x2,x3,x4), dt x3 and dx2 = x4). To do this, you will need to use Newton’s second law, which states in dt where isthevelocity(whichmeansdx1 = x2 x4 dt x1 x3 isthepositionofthepumpkinand 2 this case that where FG is the force due to gravity (which points straight down and has length m · g , where g = 9.8 m/s2) and FAR is the force due to air resistance. It’s okay to leave the constants m, g , and b as symbols.
(b) Write code to apply the (forward) Euler method to the dynamical system from part (a) in order to find out how far the pumpkin will travel if it is shot at 300 m/s at an angle of 45 degrees. Use ∆t = 0.001. Include the code, as well as how far the pumpkin went. (c) Write code which tests all angles between 20 degrees and 50 degrees (i.e. test θ = 20◦, 21◦, 22◦, ..., 50◦) and finds the angle which sends the pumpkin the farthest. Include the code as part of your answer, as well as which angle is the best and how far the pumpkin went for that angle. [Optional: plot the optimal trajectory]
Here is my code so far:
using PyPlot function simulate_orbit()
T = 1000 # total time of the simulation (in seconds)
h = 0.001 # time step size
N = int(T / h) # number of time steps (int(a) rounds a to an integer)
m = 4 # mass of pumpkin (kg)
b = 0.008 # b constant
g = 9.8 # gravity force
# create arrays to store the position and velocity at each time step
x1 = zeros(N+1, 1)
x2 = zeros(N+1, 1)
x3 = zeros(N+1, 1)
x4 = zeros(N+1, 1)
x5 = zeros(N+1, 1)
x6 = zeros(N+1, 1)
# set initial position
x1[1] = 0
x2[1] = 0
# set initial velocity
x3[1] = 300*cosd(45)
x4[1] = 300*sind(45) # Velocity of the trajectory
# forward Euler method
for n = 1:N
x1[n+1] = x1[n] + h * x3[n]
x2[n+1] = x2[n] + h * x4[n]
x3[n+1] = x3[n] + (h * (-0.6*x3[n]))
x4[n+1] = x4[n] + h * (-300*b*1/4*x4[n] - g)
end
nothing
# plot the path
x1_max = maximum(abs(x1))
x2_max = maximum(abs(x2))
a = maximum([x1_max, x2_max]) * 1.1
plot(x1, x2)
ax = gca()
ax[:set_xlim]([0, 500])
ax[:set_ylim]([0, 1000])
end
simulate_orbit()
Clearly it doesn't work and I'm not sure what I'm doing wrong. Help :( My guess is it is my x3 and x4 portion but I can't figure it out!
The force for the air resistance is quadratic in the velocity, more precisely, the accelerations compute as
Also, do a reality check. 300 m/s is fast, close to the velocity of sound. (How does the pumpkin survive the initial acceleration? It should leave the muzzle as a paste.) But not fast enough to send the pumpkin into the stratosphere. 15min simulation time is far too long, about 1min is reasonable, make it 100s to be safe.
You should somewhere in the code test for a sign change in the y coordinate and stop the simulation at that point.