# Projectile motion with air friction

I need to find the x and y components of a projectile's position, its velocity at any evaluation point, time of flight, range, and maximum height. I have already made the code for a projectile ignoring air friction, but I am having trouble figuring out how to take air friction into account. I know I need to somehow include these equations

Sum of $Fx = m \times Ax =-k\times Vx$ where $k$ is a constant

Sum of $Fy = m\times Ay =-k \times Vy-mg$

This is my code for the projectile ignoring air friction.

h = .01; g = 9.8; x = 0; y = 0; V = 40;
Vx = 40 Cos[40 °];
Vy = 40 Sin[40 °];
Vy[n_] := Vy[n - 1] - g*h;
Vx[n_] := 40 Cos[40 °];
x[n_] := x[n - 1] + 40 Cos[40 °]*h;
y[n_] := y[n - 1] + Vy[n - 1]*h;
t == (V*2*Sin[40 °])/g;
R == (V^2*Sin[2*40 °])/g;
Ymax == Sin[40 °]^2*V^2/(2 g);

• Welcome to MSE! Please consider creating an account to take advantage of the many additional features (eg. voting, favorites). – TransferOrbit Mar 22 '15 at 9:21
• I suggest you read this Wolfram Blog article – m_goldberg Mar 22 '15 at 15:47

Define constants

m = 1; g = 9.8; k = .1;


Define initial conditions

vZero = 40;
θzero = 40 Degree;
init = {
Thread[{vx, vy} == {vZero Cos[θzero], vZero Sin[θzero]}]
}


Define some force Functions

fGrav[r_List] := {0 , -m g } r
fWind[v_List] := -k v


Solve the differential equations with Solve

soln = DSolve[
Flatten@{
m {rx''[t], ry''[t]} ==
fGrav@ {rx[t], ry[t]} + fWind@ {vx[t], vy[t]}],
init
},
{rx[t], ry[t], vx[t], vy[t]},
t
];


Then plot them with ParametricPlot if you want to visualize the trajectory

ParametricPlot[Flatten@({rx[t], ry[t]} /. soln), {t, 0, 2},  AxesLabel -> {"rx", "ry"}] Remember, of course, that this is a mathematical solution. Normal ballistic objects tend not to go underground so exuberantly.

• I believe the OP wants to model a projectile launched into a tangible atmosphere. So drag as well a wind must be considered. Hitting the ground could be modeled as a discontinuity in the drag, with soil being considered a very high drag medium :-) – m_goldberg Mar 22 '15 at 18:18
• That would be a rather interesting model. I assumed here that the OP wasn't quite clear what they wanted and had just been handed a poor model of "atmospheric drag" in their problem; the fact that the OP posted a rather similar question, to which you provided a very nice answer, only decreased the weakness of my suspicion. A perusal of the excellent blog entry you posted above is well placed here. – TransferOrbit Mar 22 '15 at 19:45

I think this is a physics problem rather than a Mathematica problem. But the way you coded is problematic. You shouldn't use induction here, it will only give you integer solutions in your code. Instead, you should solve the equations.

I don't understand your question well, but I think the information I provided is sufficient.

First part, the position and velocity at any time given.

DSolve[{m*x''[t] == -k*x'[t], x == 0, x' == 40 Cos[40 °]}, x[t], t]
(*{{x[t] -> (40 E^(-((k t)/m)) (-1 + E^((k t)/m)) m Cos[40 °])/k}}*)

D[x[t] /. Flatten@%, t] // FullSimplify
(*40 E^(-((k t)/m)) Cos[40 °]*)

DSolve[{m*y''[t] == -k*y'[t] - mg, y == 0, y' == 40 Sin[40 °]}, y[t], t]
(*{{y[t] -> -((E^(-((k t)/m)) (m mg - E^((k t)/m) m mg + E^((k t)/m) k mg t +
40 k m Sin[40 °] - 40 E^((k t)/m) k m Sin[40 °]))/k^2)}}*)

D[y[t] /. Flatten@%, t] // FullSimplify
(*(-mg + E^(-((k t)/m)) (mg + 40 k Sin[40 °]))/k*)

x[t_]:= (40 E^(-((k t)/m)) (-1 + E^((k t)/m)) m Cos[40 °])/k;
Vx[t_]:= 40 E^(-((k t)/m)) Cos[40 °];
y[t_]:=(E^(-((k t)/m)) (-E^(((k t)/m))k mg t + (-1 + E^((k t)/m)) m (mg +
40 k Sin[40 °])))/k^2
Vy[t_]:=(-mg + E^(-((k t)/m)) (mg + 40 k Sin[40 °]))/k
R[t_]:=Sqrt[x[t]^2+y[t]^2]
Ymax[t_]:=(m (-mg Log[1 + (40 k Sin[40 °])/mg] + 40 k Sin[40 °]))/k^2


Second, the position and velocity at any position x=s

Solve[x[t] == s, t]
(*{{t -> ConditionalExpression[(m (2 I π C +
Log[(40 m Cos[40 °])/(-k s +
40 m Cos[40 °])]))/k, C ∈ Integers]}}*)


Simply take C=0 then plug it into y[t], you will get

y[s_]:=(mg (-40 m Log[1 + (k s)/(-k s + 40 m Cos[40 °])] +
k s Sec[40 °]))/(40 k^2) + s Tan[40 °]


for y=h, you can follow the similar steps.