Prev: Solving a transcendental equation
Next: unit_impulse, unit_step and unit_ramp functions
From: Brian Ickler on 31 Jan 2010 18:36
Hello -
Using ode45, I can plot the trajectory of a projectile using basic equations of motion (no drag) as follows (I've spared you some of the initial calculations and plot formatting):

global g v0x v0y
g = 32.2; % gravity acceleration
v0x = 103.9; % initial x velocity
v0y = 60.0; % initial y velocity
x0 = 0; % initial x position
y0 = 0; % initial y position

tspan = [0, 5];
init_cond = [x0; y0];
options = odeset('Events', @events);
[t,y,te,ye,ie] = ode45(@rk4_nodrag, tspan, init_cond, options);
xdir = y(:,1);
ydir = y(:,2);
plot(xdir,ydir)

% dydt(1) is dx/dt = x direction velocity, dydt(2) is dydt = y direction velocity
function dydt = rk4_nodrag(t,y)
global g v0x v0y
dydt = [v0x
v0y - (g * t)];
end

The problem is when I change the system of equations to include drag force, which includes a velocity^2 term in each direction. As I have it coded below, I am using v0x and v0y, whereas the drag force should be calculated at each time step, using the current velocity, not the initial velocity.

% dydt(1) is dx/dt = x direction velocity, dydt(2) is dydt = y direction velocity
function dydt = rk4_drag(t,y)
global g v0x v0y m Cd rho A
dydt = [v0x - (1/(2 * m) * Cd * rho * A * v0x^2 * t)
v0y - (g * t) - ((1/(2 * m)) * Cd * rho * A * v0y^2 * t)];
end

Any help? I've tried, but failed to change v0x^2 and v0y^2 to something else that would be the velocity (dydt(1) and dydt(2)) at the current position of the projectile.

Many thanks,
Brian
 | 
Pages: 1
Prev: Solving a transcendental equation
Next: unit_impulse, unit_step and unit_ramp functions