# Projectile motion in Mathematica using Euler's method [duplicate]

I need to write the equations of a ball in projectile motion (ignoring air friction) with an initial velocity of 40 m/s at an angle of 40° with respect to the horizontal. Specifically, I need to write:

1. x (n) and y (n) for the projectile (where n represents the nth evaluation point)
2. The x and y components of velocity (Vx(n) and Vy(n))
3. the time of flight, range and maximum height achieved (and compare to the theoretical values) (I'm not sure how to write these at all)

This is what I have in my program so far, however when I compile testing Vy[1], I'm not sure what the answers at the bottom are referring to.

h = .01; g = 9.8; x[0] = 0; y[0] = 0;
Vy[0] = 40 Sin[40 \[Degree]]
Vy[n_] := Vy[n - 1] - g*h
Vx[n_] := 40 Cos[40 \[Degree]]
x[n_] := x[n - 1] + 40 Cos[40 \[Degree]]*h
y[n_] := y[n - 1] + Vy[n - 1]*h
Vy[1]

• Note that Mathematica's trig functions are defined in terms or radians, so to use 40° as an argument, you must write Sin[40 Degree] or Sin[40 °]. – m_goldberg Mar 21 '15 at 23:11

I will leave to you work out the theoretical solution and just present a way to make the estimates of the flight time, the maximum range, and the maximum height.

I will start with a rather crude implementation of Euler's method implemented as a simple While-loop.

pts = {};
Module[{t = 0., dt = .01, a = -9.807, x = 0., y = 0.,
vx = 40. Cos[40. Degree], vy = 40. Sin[40. Degree]},
While[y >= 0.,
t += dt;
vy += a dt; (* don't need to update vx -- it is constant *)
x += vx dt; y += vy dt;
pts = {pts, {t, x, y}}];
pts = Partition[Flatten[pts], 3]];


The list pts accumulates the trajectory points in the form {t, x, y}. To examine both ends of the list, I will use Short.

Short[pts, 3]

{{0.01, 0.306418, 0.256134}, <<522>>, {5.24, 160.563, -0.167002}}


From the above it already clear that the ball hit the ground about 160.5 meters from where it was thrown in about 5.25 seconds.

To get a better feel for the results, I plot them (I always do this after running a model).

ListLinePlot[{pts[[All, ;; 2]], pts[[All, 1 ;; 3 ;; 2]]},
AxesLabel -> {"t", "range / height"},
PlotLegends -> LineLegend[{"range: x[t]", "height: y[t]"}]]


The plot looks good -- that is, it looks as I expected it to, so I can usingproceed to finding the flight time, the maximum range, and the maximum height. The easiest way (in the sense of having Mathematica do all the drudge work) to get these values form the list of points is to derive interpolating functions from them.

range = Interpolation[pts[[All, ;; 2]]]; (* interpolates x[t] *)
height = Interpolation[pts[[All, 1 ;; 3 ;; 2]]]; (* interpolates y[t] *)


Now I solve for the required values.

Off[InterpolatingFunction::dmval]; (* don't want warning *)
flightTime = t /. FindRoot[height[t] == 0., {t, 5.}]

5.2335 (* seconds *)

maxRange = range[flightTime]

 160.364 (* meters *)


As the plot shows, the height curve is symmetric about the middle of the flight time and, thus, reaches its maximum there.

maxHeight = height[flightTime/2]

 33.5761 (* meters *)


Please be more specific, what do you mean by Euler method, because normally it refers to Finite Difference Method for Differential Equation. So do you really intend to solve m*x'' = F? Or do you just want to use the standard approach?:

1. Determine the time of flight by in y direction (up)
2. Applying the fact that there is a constant motion in the x direction

Next, @m_goldberg mentioned, you are using wrong units. NEVER use subscripts and/or superscripts in Mathematica for the name of variables, functions etc. Use them only if builtin symbols require that (Sum). Mathematica is unable to process that and strange results may happen and this is the reason for the first weird result (40Sin[40]), not the fact that you are not using radians as the Sin's domain is real numbers. Further, I am not sure what you mean by

x[n_] := x[n] = ...
y[n_] := y[n] = ...


Here I made minor changes to your code to produce what you ask - v[1] according to your definition:

h = 0.01;
g = 9.8;
x[0] = 0;
y[0] = 0;
v0 = 40;
\[Alpha]0 = 40*\[Pi]/180;
vy[0] = v0*Sin[\[Alpha]0];
vy[n_] := vy[n - 1] - g*h;
x[n_] := x[n - 1] + v0*Cos[\[Alpha]0]*h;
y[n_] := y[n - 1] + vy[n - 1]*h;
vy[1]
Out[84]= 25.6135


Please try to specify what you are trying to actually accomplish because it is not clear from the code.