-1
$\begingroup$

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] = 0; y[0] = 0; V[0] = 40;
Vx[0] = 40 Cos[40 °];
Vy[0] = 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[0]*2*Sin[40 °])/g;
R == (V[0]^2*Sin[2*40 °])/g;
Ymax == Sin[40 °]^2*V[0]^2/(2 g);
| improve this question | | | | |
$\endgroup$
  • $\begingroup$ Welcome to MSE! Please consider creating an account to take advantage of the many additional features (eg. voting, favorites). $\endgroup$ – TransferOrbit Mar 22 '15 at 9:21
  • 1
    $\begingroup$ I suggest you read this Wolfram Blog article $\endgroup$ – m_goldberg Mar 22 '15 at 15:47
1
$\begingroup$

Define constants

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

Define initial conditions

vZero = 40;
θzero = 40 Degree;
init = {
  Thread[{rx[0], ry[0]} == {0, 0}],
  Thread[{vx[0], vy[0]} == {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@{
     Thread[
      m {rx''[t], ry''[t]} == 
       fGrav@ {rx[t], ry[t]} + fWind@ {vx[t], vy[t]}],
     Thread[{vx[t], vy[t]} == {rx'[t], ry'[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"}]

Blockquote

Remember, of course, that this is a mathematical solution. Normal ballistic objects tend not to go underground so exuberantly.

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ 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 :-) $\endgroup$ – m_goldberg Mar 22 '15 at 18:18
  • $\begingroup$ 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. $\endgroup$ – TransferOrbit Mar 22 '15 at 19:45
1
$\begingroup$

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] == 0, x'[0] == 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] == 0, y'[0] == 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[1] + 
   Log[(40 m Cos[40 °])/(-k s + 
     40 m Cos[40 °])]))/k, C[1] ∈ Integers]}}*)

Simply take C[1]=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.

| 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.