2
$\begingroup$

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]
$\endgroup$
  • 1
    $\begingroup$ Try this ParametricPlot[{x1[t], y1[t]} /. Movement2, {t, 0, 2}] $\endgroup$ – zhk May 2 '17 at 17:44
  • $\begingroup$ Isn't final time missing? $\endgroup$ – zhk May 2 '17 at 17:53
  • $\begingroup$ 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 :) $\endgroup$ – Andrej Licanin May 2 '17 at 18:15
  • $\begingroup$ 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 . $\endgroup$ – Andrej Licanin May 2 '17 at 18:25
  • $\begingroup$ Two different values for k. Why you expect that the two output should be the same ? $\endgroup$ – zhk May 2 '17 at 18:35
2
$\begingroup$

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]

enter image description here

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?

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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 :) $\endgroup$ – Andrej Licanin May 2 '17 at 20:12
  • $\begingroup$ @AndrejLicanin Please, see my edit. $\endgroup$ – yarchik May 2 '17 at 20:22
0
$\begingroup$

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}]   

enter image description here

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

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.