Projectile trajectory with quadratic air resistance in MATLAB

[Hoppa]

this is the problem:


A projectile is fired with an initial speed u and at an angle of elevation ®. The air
resistance is known to be quadratic and the terminal velocity has a magnitude vt. Show
that the equations of motion for the projectile can be cast into the form
y' = f (t; y); t >= 0

Using MATLAB solve the differential equations for an initial speed of 200 ms¡1,
an angle of elevation of 45± and a terminal velocity of 250 ms¡1. Plot the trajectory of
the projectile and, on the same graph, plot the trajectory that the projectile would have
in the absence of air resistance. Use a time range that allows both trajectories to at least
return to their initial height, without going significantly beyond that position.

 

[Hoppa]

and here is what I've got so far..

Initial speed = m
Angle of elevation = a

y = é x ù é vx ù é 0 ù
| vx ô , f(t,y) = | -gvt-2vx Övx2 + vz2 | , y(0) = | m cos a |
| z | | vz | | 0 |
ë vz û ë -g(1 + vt-2vzÖvx2 + vz2 û ë m sin a û

Where x and z are the horizontal and vertical coordinates and vx and vz are the corresponding velocities.

The position of the projectile at time t is given by:

x = (m cos a ) t ( 1 – e–t/t )

y = -gvtt + t (m sin a + gvt) (1 – e–t/t )

Where gvt is the magnitute of the terminal velocity.


Matlab code

global vterm tau v0y
g = 9.8;
vterm = 250;
tau = vterm/g;
v0 = 200;
alpha = 45;
time = 0:2:50
range = time;
n = length(time);
for i = 1:n
angle = alpha(i)*pi/180;
v0x = v0*cos(angle);
v0y = v0*sin(angle);
time = fzero(’yproject’,[10;50]);
range(i) = v0x*tau*(1-exp(-time/tau));
end
plot(alpha, range);
xlabel(’Elevation (degrees)’);
ylabel(’Range (metres)’);
which also uses the MATLAB function:
function ypos = yproject(t)
global vterm tau v0y
ypos = -vterm*t + tau*(v0y+vterm)*(1-exp(-t/tau));

 

[thepatientmental]

To solve this problem, we need to first understand the equations of motion for a projectile with quadratic air resistance. These equations can be derived from the basic principles of motion and are given by:

x' = u cos(θ) - (k/m)vx|vx|
y' = u sin(θ) - g - (k/m)vy|vy|

where x and y are the position coordinates, u is the initial speed, θ is the angle of elevation, k is a constant related to air resistance, m is the mass of the projectile, vx and vy are the velocities in the x and y directions respectively, and g is the acceleration due to gravity.

We can rewrite these equations in terms of the time derivative of the position coordinates, denoted by x' and y', which gives us:

x'' = - (k/m)vx|vx|
y'' = - g - (k/m)vy|vy|

Now, we can use MATLAB to solve these differential equations numerically. We start by defining the variables and constants:

u = 200; % initial speed in m/s
θ = 45*pi/180; % angle of elevation in radians
vt = 250; % terminal velocity in m/s
k = 0.01; % air resistance constant
m = 1; % mass of the projectile in kg
g = 9.8; % acceleration due to gravity in m/s^2

Next, we define the differential equations in terms of MATLAB's anonymous function notation:

f = @(t,y) [y(3); y(4); -(k/m)*y(3)*abs(y(3)); -g-(k/m)*y(4)*abs(y(4))];

Here, y(1) and y(2) represent the x and y coordinates respectively, while y(3) and y(4) represent the velocities in the x and y directions respectively.

Finally, we can use MATLAB's built-in ode45 function to solve the differential equations and plot the trajectories:

[t,y] = ode45(f, [0 10], [0 0 u*cos(θ) u*sin(θ)]); % time range is chosen as 0 to 10 seconds
plot(y(:,1), y(:,2)); % plot the trajectory with air resistance
hold on
plot(u*cos(θ)*t, u*sin(θ)*t - 0.5*g