-3
$\begingroup$

This question already has an answer here:

I am trying to find the time and maximum height of a projectile motion with air friction however Mathematica is "running" for a really long time and not giving an answer. Is there something wrong with my code?

Clear[X, Vx, Y, Vy]
V[0] = 40; Y[0] = 0; h = .01; m = 1; k = .5; g = 9.8; Ax[0] = 0; 
Vx[0] = 40 Cos[40 \[Degree]]; Vy[0] = 40 Sin[40 \[Degree]]
X[n_] := X[n] = X[n - 1] + V[n - 1]*h
Vx[n_] := Vx[n] = Vx[n - 1] + Ax[n - 1]*h
Ax[n_] := (-k*Vx[n - 1])/m
Y[n_] := Y[n] = Y[n - 1] + Vy[n - 1]*h
Vy[n_] := Vy[n] = Vy[n - 1] + Ay[n - 1]*h
Ay[n_] := ((-k*Vy[n - 1])/m) - g
t = Catch[Do[If[Y[n] < 0, Throw[n - 1]], {n, 1000}]]
Ymax = Catch[Do[If[Vy[n] < 0, Throw[Y[n - 1]]], {n, 1000}]]
$\endgroup$

marked as duplicate by Kuba, bbgodfrey, Sjoerd C. de Vries, Dr. belisarius, Jens Mar 26 '15 at 0:42

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • $\begingroup$ You haven't defined Ay[0] and V[n_] yet. $\endgroup$ – Taiki Mar 23 '15 at 16:35
  • $\begingroup$ It runs for a short time and give hints (messages) as to what is wrong. What version are you running? It's hard to imagine any version in which $RecursionLimit behaves differently than in V10.0.2. $\endgroup$ – Michael E2 Mar 23 '15 at 23:33
  • $\begingroup$ I am using version 9.0. I have waited 10+ minutes and still am not getting an error message, even after defining Ay[0] $\endgroup$ – user700 Mar 24 '15 at 3:05
  • $\begingroup$ This is a thoroughly glorious use of Catch and Throw. However, I think a simpler algorithm, possible using Do or Table or Scan might be what you need. $\endgroup$ – djp Mar 24 '15 at 9:48
  • $\begingroup$ For its calculation Vx[1] needs Ax[0], which, on its turn, needs Vx[-1] which is undefined $\endgroup$ – Sjoerd C. de Vries Mar 25 '15 at 19:46
3
$\begingroup$

There is something wrong with your code-- it has terrible "big-O" complexity. Every time you call Y[n], it has to go and recalculate every single preceding value of Y[n-1] etc., as well as all the other variables.

First@AbsoluteTiming[Y[29]]
(* 16.3 seconds *)
First@AbsoluteTiming[Y[30]]
(* 26.48 seconds *)

And you're trying to calculate all values of Y potentially up to 1000...

If you use memoization it will sing along. Make the following changes:

X[n_] := X[n] = X[n - 1] + V[n - 1]*h
Vx[n_] := Vx[n] = Vx[n - 1] + Ax[n - 1]*h
Y[n_] := Y[n] = Y[n - 1] + Vy[n - 1]*h
Vy[n_] := Vy[n] = Vy[n - 1] + Ay[n - 1]*h

First@AbsoluteTiming@Y[1000]
(* 0.017 seconds - the first time.  *)

Better be careful with memoization, because it's easy to re-use incorrect old definitions. You should use Clear[X, Vx, Y, Vy] regularly.

Throw...Catch ... So many cleaner ways to do it, but I can't think of any so amusing.

$\endgroup$
  • $\begingroup$ I tried your changes but still am getting an error...I also tried other ways besides catch and throw but it still didnt work $\endgroup$ – user700 Mar 24 '15 at 22:44
  • $\begingroup$ Mate you need to refine your question. What is your question about Throw/Catch? If your question is not about Throw/Catch, you need to edit it and make it specific. I'm happy to help you debug code if you reciprocate. $\endgroup$ – djp Mar 25 '15 at 6:33
  • $\begingroup$ You reverted to the broken version where Ay[0] wasn't defined? $\endgroup$ – djp Mar 25 '15 at 7:43
2
$\begingroup$

Or use RSolve

h = 1/100; m = 1; k = 1/2; g = 98/10; V0 = 40;

sys = RSolve[{
     X[n] == X[n - 1] + Vx[n - 1]*h,
     Vx[n] == Vx[n - 1] + Ax[n - 1]*h,
     Ax[n] == (-k*Vx[n - 1])/m,
     Y[n] == Y[n - 1] + Vy[n - 1]*h,
     Vy[n] == Vy[n - 1] + Ay[n - 1]*h,
     Ay[n] == ((-k*Vy[n - 1])/m) - g,
     X[0] == 0,
     Y[0] == 0,
     Vx[0] == V0 Cos[40 Degree],
     Vy[0] == V0 Sin[40 Degree],
     Ax[0] == 0,
     Ay[0] == g},
    {X[n], Vx[n], Ax[n], Y[n], Vy[n], Ay[n]},
    n][[1]] // FullSimplify

{Ax[n] -> 100/7 Sqrt[ 2] ((1/20 (10 - 7 Sqrt[2]))^n - (1/20 (10 + 7 Sqrt[2]))^n) Cos[ 40 [Degree]], Vx[n] -> 1/7 2^(3/2 - 2 n) 5^( 1 - n) (-(10 - 7 Sqrt[2])^(1 + n) + (10 + 7 Sqrt[2])^(1 + n)) Cos[ 40 [Degree]], X[n] -> 1/7 2^(1 - 2 n) 5^-n (7 2^(3 + 2 n) 5^(1 + n) + (10 - 7 Sqrt[2])^ n (-140 + 99 Sqrt[2]) - (10 + 7 Sqrt[2])^n (140 + 99 Sqrt[2])) Cos[ 40 [Degree]], Ay[n] -> 1/7 2^(-(3/2) - 2 n) 5^(-1 - n) ((10 + 7 Sqrt[2])^ n (49 (-30 + 7 Sqrt[2]) - 2000 Sin[40 [Degree]]) + (10 - 7 Sqrt[2])^ n (49 (30 + 7 Sqrt[2]) + 2000 Sin[40 [Degree]])), Vy[n] -> 1/7 2^(-(3/2) - 2 n) 5^(-2 - n) (-343 2^(5/2 + 2 n) 5^(1 + n) + (10 - 7 Sqrt[2])^ n (49 (-101 + 70 Sqrt[2]) + 1000 (-10 + 7 Sqrt[2]) Sin[40 [Degree]]) + (10 + 7 Sqrt[2])^ n (49 (101 + 70 Sqrt[2]) + 1000 (10 + 7 Sqrt[2]) Sin[40 [Degree]])), Y[n] -> 1/7 4^(-1 - n) 5^(-3 - n) (-343 2^(1 + 2 n) 5^n n + 7 2^(1 + 2 n) 5^n (9849 + 20000 Sin[40 [Degree]]) + (10 - 7 Sqrt[2])^ n (49 (-1407 + 995 Sqrt[2]) + 1000 (-140 + 99 Sqrt[2]) Sin[40 [Degree]]) - (10 + 7 Sqrt[2])^ n (49 (1407 + 995 Sqrt[2]) + 1000 (140 + 99 Sqrt[2]) Sin[40 [Degree]]))}

Plot[Evaluate[{X[n], Y[n]} /. sys],
 {n, 0, 403}, PlotLegends -> {X, Y}]

enter image description here

$\endgroup$

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