# Projectile motion with drag

first post here :)

So i have a problem to solve projectile motion with drag when drag is linear and quadratic

This is my first take

m x''[t] == -k x'[t]
m y''[t] == -m g - k y'[t]
Movement =
DSolve[{m x''[t] == -k x'[t], m y''[t] == -m g - k y'[t], x[0] == 0,
y[0] == 0, x'[0] == Vo Cos[\[Alpha]],
y'[0] == Vo Sin[\[Alpha]]}, {x[t], y[t]}, t]

Data = {m -> 1, Vo -> 20, g -> 9.81, k -> 0.5, \[Alpha] -> \[Pi]/6,
tt -> 2}
ParametricPlot[{x[t] /.Movement, y[t] /. Movement} /. Data, {t, 0,
tt /. Data}]


It gives a nice graph but im not sure if i made the right assumption for drag-that i can just put it like x'[t] and y[t] in the x and y axis .

This is my second take with

v=sqrt[x1'[t]^2 +y1'[t]^2]
m x1''[t] == -k Sqrt[x1'[t]^2 + y1'[t]^2]
m y1''[t] == -m g - k Sqrt[x1'[t]^2 + y1'[t]^2]
Data2 = {m -> .1, Vo -> 10, g -> 9.81,
k -> 0.001, \[Alpha] -> \[Pi]/6, ttt -> 2}
Movement2 =
NDSolve[{m x1''[t] == -k Sqrt[x1'[t]^2 + y1'[t]^2],
m y1''[t] == -m g - k Sqrt[x1'[t]^2 + y1'[t]^2], x1[0] == 0,
y1[0] == 0, x1'[0] == Vo Cos[\[Alpha]],
y1'[0] == Vo Sin[\[Alpha]]} /. Data2, {x1, y1}, {t, 0, 100}]
Plot[{x1[t] /. Movement2, y1[t] /. Movement2}, {t, 0, 2}]


This gives me a solution but when i plot it i dont know what am i reading. How do i combine the solution from NDSolve (x and y ) to get a single line trajectory.(tried parametric plot,didnt work) Also witch one of these 2 takes is right ?

Edit:Fixed and found the solution for linear drag,now when i enter for quadratig im getting trouble with singularity in NDSolve .

Data = {m -> 1, Vo -> 20, g -> 9.81, k -> 0.5, \[Alpha] -> \[Pi]/6,
tt -> 2}
Movement3 =
NDSolve[{m x1''[t] == -k x1'[t]^2, m y1''[t] == -m g - k y1'[t]^2,
x1[0] == 0, y1[0] == 0, x1'[0] == Vo Cos[\[Alpha]],
y1'[0] == Vo Sin[\[Alpha]]} /. Data, {x1, y1}, {t, 0, 100},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4},
WorkingPrecision -> 50, AccuracyGoal -> 15, PrecisionGoal -> 15,
MaxSteps -> 20000]
gr2 = ParametricPlot[{x1[t], y1[t]} /. Movement2, {t, 0, 1}]


I keep getting a singularity that changes as i change the paramateres.Can i fix this somehow? At around 1.217 it hits the ground but when i go beyond that (lets say {t,0,2}) the graph goes bonkers.

And when i try DSolve i get a nicer solution but i have imaginary numbers,how do i plot this solution?

Movement3 =
DSolve[{m x1''[t] == -k x1'[t]^2, m y1''[t] == -m g - k y1'[t]^2,
x1[0] == 0, y1[0] == 0, x1'[0] == Vo Cos[\[Alpha]],
y1'[0] == Vo Sin[\[Alpha]]}, {y1[t], x1[t]}, t]

• Try this ParametricPlot[{x1[t], y1[t]} /. Movement2, {t, 0, 2}]
– zhk
Commented May 2, 2017 at 17:44
• Isn't final time missing?
– zhk
Commented May 2, 2017 at 17:53
• Lol it worked,i used it with /.Movement2 inside curly's next to x1 and y1 before. Final time u mean time of landing? havent gotten around to that i wanna get the differential eqations right for both cases. Thank you :) Commented May 2, 2017 at 18:15
• Well,i managed to graph the second one,and its not the same as first one there is a notable difference on range of the projectile,with same data .Still dont know witch one is right . Commented May 2, 2017 at 18:25
• Two different values for k. Why you expect that the two output should be the same ?
– zhk
Commented May 2, 2017 at 18:35

Your second equation is not completely correct. The declaration is proportional to the velocity squared, i.e., $a_d=-kv^2$ but this needs to be properly projected. This means $\boldsymbol{a}_d=-k\,a_d\, \boldsymbol{v}/v=-k\boldsymbol{v} v$ because drag is always oposite to the velocity direction, which is $\boldsymbol{v}/v$.

v = Sqrt[x'[t]^2 + y'[t]^2];
ics = {x[0] == 0, y[0] == 0, x'[0] == Vo Cos[α],
y'[0] == Vo Sin[α]};
eqs[1] = {m x''[t] == -k x'[t], m y''[t] == -m g - k  y'[t]};
eqs[2] = {m x''[t] == -k v x'[t], m y''[t] == -m g - k v y'[t]};
eqs[3] = {m x''[t] == -k v , m y''[t] == -m g - k v };
sys[i_] := Join[eqs[i], ics, {WhenEvent[y[t] <= 0, y'[t] -> 0.01]}];
data = {m -> .1, Vo -> 100, g -> 9.81,
k -> 0.001, α -> Pi/6};
move = Table[
First@NDSolve[sys[i] /. data, {x, y}, {t, 0, 10}], {i, 3}];
gr = ParametricPlot[Evaluate[({x[t], y[t]} /. move)], {t, 0, 10},
AspectRatio -> 1, PlotLegends -> Automatic]


The 1st case is drag linearly proportional to the velocity, 2nd is the quadratic drag, while 3rd is your original formula. Finally, the particle collides with the "ground". Therefore the WhenEvent.

I wonder though why AspectRatio is not working as expected and why the propagation slows down if in WhenEven I reduce the reflection velocity to 0.00001. Can someone improve my answer?

• Thank you,so when you project the quadratic drag on x it goes -k * [v^2 * x'[t])^2 * x'[t]/v ) ,the x'[t]/v is the cos of angle for x axis ,?just checking need to be sure :) Commented May 2, 2017 at 20:12
• @AndrejLicanin Please, see my edit. Commented May 2, 2017 at 20:22

If you look at your Movement2 variable, it is a 1x2 array. Then you look at x1[t] /. Movement2, it returns a list. Plot[ ] handles this, ParametricPlot[ ] does not. So just grab the first (and only) item from the list. Add this to your code above:

ParametricPlot[{(x1[t] /. Movement2) // First, (y1[t] /. Movement2) //
First}, {t, 0, 2}]


EDIT: although I like @zhk 's answer better. In general, got to keep an eye on your list-ology.