# Projectile motion with air resistance proportional to velocity squared, system of DE's.

A plane flies at a constant altitude of 1000 ft with a constant speed of 300 mph. The plane drops a relief package to a person on the ground. Assume the origin is point where the supply pack is released and the positive x-axis points forward from the plane, while the positive y-axis points downward. Assume the horizontal and vertical components of air resistance are proportional to the square of the velocity. Assume that the constant of proportionality for the air resistance is k = 0.0053 and that the package weighs 256 lbs.

Find the horizontal distance the package will travel from the time of its release to the point where it hits the ground.

I realize I need to build a system of differential equations and will need to solve this numerically, but I am not sure I am building my equations correctly.

I think I need to consider the following DE for each direction.

$m\frac{d^2x}{dt^2}=mg-k(\frac{dx}{dt})^{2}$

Given that

$x(0)=0$, $x'(0)=300$, $y(0)=0$, $y(0)=0$, $m=\frac{w}{g}=\frac{256}{32}=8$

I built the equations as follows:

$8x''= -.0053(x')^2$

$8y''=256-.0053(y')^2$

I cannot get a solution from wxMaxima, therefore I assume I am building the equations incorrectly.

Can I get some help on setup with this system of DE's?

• Why would there be a gravitational force in the horizontal direction? I am also a bit confused by your weight/mass units... – Ian Apr 22 '18 at 10:02
• Ahh, there wouldn't. So the equation for horizontal direction should be $8x'' = -0.0053(x')^{2}$? – PerpetualStudent Apr 22 '18 at 10:07
• The weight of the package is 256lbs and so if $m=\frac{W}{g}$ where $g = 32ft/sec^{2}$, then $m=\frac{256}{32}=8 slug$, correct? – PerpetualStudent Apr 22 '18 at 11:19
• @Denise : No, the air resistance couples both equations, the force is $-k\|(\dot x,\dot y)\|_2\,(\dot x,\dot y)$. – Dr. Lutz Lehmann Apr 22 '18 at 13:22
• Ok, it is not easy to find sources about air resistance in components, as that is only numerically solvable and does not make good manual examples. Per the text, the resistance is proportional to the square of the velocity, $|F_{air}|=k\,|v|^2$. Now the direction is opposite to the velocity, which is $(-v_x/|v|,\,-v_y/|v|)$. Together that makes the directed force $\vec F_{air}=(-k\,|v|\,v_x,\,-k\,|v|\,v_y)$. – Dr. Lutz Lehmann Apr 22 '18 at 13:57

A better model could be done

$$m \ddot x = - k \sqrt{\dot x^2+ \dot y^2}\dot x\\ m \ddot y = -m g - k \sqrt{\dot x^2+ \dot y^2}\dot y\\$$

Integrating with parameters $m = 1, k = 0.005, v_{x_0} = 10, v_{y_0} = 20, x_0 = y_0 = 0, g = 10$ the following plot can be obtained

In red without the aerodynamic force and in black with the aerodynamic force

• Yes I agree with you. My fault. Your model is more accurate. I will correct my answer accordingly. – Cesareo Apr 22 '18 at 13:32

You need to first convert everything to use the same units, here feet and pound and seconds. Then $x(0)=0$, $y(0)=1000ft$, $x'(0)=300\,\frac{5280ft}{3600s}=440\frac{ft}{s}$, $y'(0)=0$, $g=32\frac{ft}{s^2}$ and presumably $k= 0.0053\frac{lbs}{ft}$.

g, k, m = 32, 0.0053, 256
def odesys(u,t):
x,y,vx,vy = u
v = np.hypot(vx,vy)
return [ vx, vy, -k/m*v*vx, -g-k/m*v*vy ]
def report(u,t):
print "t=%8.6f (x,y)=(%12.10f, %12.10f)%(t,u[0],u[1])

u = [ 0, 1000, 440, 0 ]
dt = 1;
t = 0;
while u[1] > 0:
u1 = odeint(odesys, u, [t, t+dt])[-1]
if abs(u[1]) < 1e-6: break
if u1[1]>0:
u, t = u1, t+dt
report(u,t)
else:
dt = -u[1]*2/(u[3]+u1[3])


with the resulting table

t=1.000000 (x,y)=(438.0071507243, 984.0482925753)
t=2.000000 (x,y)=(872.0659485475, 936.3846433523)
t=3.000000 (x,y)=(1302.2161300210, 857.2945582413)
t=4.000000 (x,y)=(1728.4769062757, 747.0648211872)
t=5.000000 (x,y)=(2150.8486364701, 605.9875302307)
t=6.000000 (x,y)=(2569.3147184315, 434.3637374476)
t=7.000000 (x,y)=(2983.8435931887, 232.5064350959)
t=8.000000 (x,y)=(3394.3908256516, 0.7430872659)
t=8.002844 (x,y)=(3395.5526359051, 0.0416630833)
t=8.003013 (x,y)=(3395.6216203863, 0.0000066767)
t=8.003013 (x,y)=(3395.6216314429, 0.0000000001)