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

Question

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?

Answer 1

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

Answer 2

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 speed $|v|=\sqrt{x'^2+y'^2}$, $|F_{\rm air}|=k|v|^2$. Now the direction is opposite to the velocity vector $v=(x',y')$, which is $(−x'/|v|,−y'/|v|)$. Together that makes the directed force $$ F_{\rm air}=(−k|v|x',−k|v|y'). $$ To start with the computation, 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])

where in the last step with the zero crossing a Newton-like formula using $y(t)$ and $y'(t+h/2)\approx (y'(t)+y'(t+h/2)$ is used to get a better value for the event time, then backtracking and integrating to that point gives corrected values for that estimate, repeat until the error tolerance is met.

This results in the 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)