Lab 2 - Motion In 2D

1. Write down each algebraic equation you used calculate these results. First write an algebraic expression for \(t\), then write one for \(v_{ix}\). You will use these equations to propagate error. \[t(y) = \qquad\qquad\qquad v_{ix}(x_f,t)=\]

2. Use the error propagation formula to find the propagated error in time

   \[E[t] = \sqrt{\left( \frac{\partial t}{\partial y} \cdot E[y] \right)^2 } = \frac{\partial t}{\partial y} \cdot E[y]\]

3. Use the error propagation formula to find the propagated error in velocity

   \[E[v_{ix}] = \sqrt{\left( \frac{\partial v_{ix}}{\partial t} \cdot E[t] \right)^2 + \left( \frac{\partial v_{ix}}{\partial x_f} \cdot E[x_f] \right)^2 }\]

Checkpoint:

• Complete the derivation of error formulae for time and velocity
• Propagate error in velocity \(E[v_{ix}]\)

Propagating error leads to complications even in simple cases, such as kinematics. If you attempt to go headlong into the math, you arrive at the following: (This includes a clever substitute-eliminate trick)

Process to solve for range in 2D free fall without knowing time:

• write out the kinematic equations \[x_f = x_i + v_i cos(\theta)\,t\] \[y_f = y_i + v_i sin(\theta)\,t + \frac12 a \, t^2\]
• solve for time with the x equation \[t=\frac{x_f - x_i}{v_{i}cos(\theta)}\]
• substitute this result into \(y_f\) equation \[y_f = y_i + v_i sin(\theta)\left(\frac{x_f - x_i}{v_{i}cos(\theta)}\right) + \frac12 a \, \left(\frac{x_f - x_i}{v_{i}cos(\theta)}\right)^2\] \[y_f = y_i + (x_f - x_i) tan(\theta) + \frac12 a \, \left(\frac{x_f - x_i}{v_{i}cos(\theta)}\right)^2\]

this is a useful result to use for solving physics problems (this technique is used for the first exam). But realize the part we are after for this experiment is actually the \((x_f - x_i\)) term, which is quadratic. Yikes! This is a show stopper and other methods must be used (sampling or other simulations).

The ’simplest’ sampling method is usually a Monte Carlo type simulation. Monte Carlo simply uses raw compute power and random numbers to brute force a result. In our case, this means inputting our initial conditions and running a script which calculates the posterior distribution of landing points. Lets start with a simple example of monte carlo: Calculating the number \(\pi\):

Lets now use the Monte Carlo method to simulate our ball launch.

Simply using the kinematics equations and giving normally distributed random values for the initial height, initial velocity, and angle of incline.

from numpy.random import normal # import Gaussian RNG

# Declare measured initial conditions
yi, Eyi = 1, 0.002
vi, Evi = 4.6, 0.1
theta, Etheta = 30*numpy.pi/180, 1*numpy.pi/180
ay = -9.8027

# create a list of uniforly spaced numbers for the time axis
time=numpy.array(numpy.arange(0,1,0.01))

# perform the kinematics calculation a bunch o times
for i in range(100):
  # x=xi+vix*t = 0 + vi*sin(theta)*t
  xpath=normal(vi,Evi)*numpy.cos(normal(theta,Etheta))*time
  # y=yi+viy*t-4.9t^2 = yi + vi*sin(theta)*t-49.t^2
  ypath=normal(yi,Eyi)+normal(vi,Evi)*numpy.sin(normal(theta,Etheta))*time+0.5*ay*time**2

  plt.plot(xpath,ypath)

plt.ylim(0,1.4)
plt.xlabel('Horizontal Position (m)')
plt.ylabel('Vertical Position (m)')
plt.title('Projectile Motion')
plt.gca().set_aspect('equal')

plt.savefig("img/02-projectile.png",bbox_inches='tight')
plt.close()
#plt.show()

This is a very illustrative figure, but the code is not calculating what we want (you need an estimate of the horizontal range to report to the boss in an ’abstract’). The above code calculates x and y for a list of times (perfect for making a plot of the trajectories) but we want a distribution of horizontal ranges at a specific value of height. To proceed with the above data, you could search thru the list and find the list index (corresponding to a time) which the vertical position becomes negative. This sounds like a fine exercise for a programming assignment, but I wish to proceed in a different direction.

Instead, follow the standard method of calculating the range of an inclined shot from physics theory:

1. Break up the velocity vector into components
2. Use the y direction to solve for time (quadratic equation)
3. Use the calculated time to find the range

Lets write a function to evaluate the solutions to the kinematic equations:

from numpy.random import normal # import Gaussian RNG
# vars to use as input: vi, theta, yi

# Step 1: vix=vi*cos(theta) viy=vi*sin(theta)
# Step 2: Kinematics in y-direction: 0=yi+vi*sin(theta)*t - 4.9 t^2
# this is a qadratic equation with a=-4.9 b=vi*sin(theta) c=yi
# t = (-b - sqrt(b*b - 4*a*c))/(2*a)
# Step 3: Kinematics x-direction: x = xi + vi*cos(q)*t

# requires 3 inputs: vi, theta, yi
# returns the calculated horizontal range
def ballistic_sim(vi,theta,yi):
  a,b,c=-4.9, vi*numpy.sin(theta), yi
  t = (-b - numpy.sqrt(b*b - 4*a*c))/(2*a)
  return vi*numpy.cos(theta)*t

# Calculate the trajectory for a randomly chosen set of input parameters
# uses numpy.random.normal() to generate random input values for vi, theta, and yi
for i in range(10):
  print("random shot",i,":",ballistic_sim(normal (vi,Evi),normal(theta,Etheta),normal(yi,Eyi)) )

random shot 0 : 2.8253864893107106
random shot 1 : 2.943230448713085
random shot 2 : 2.9486365004658768
random shot 3 : 2.942945481287842
random shot 4 : 3.0416240128377066
random shot 5 : 2.9694014633532553
random shot 6 : 2.889390480484175
random shot 7 : 2.8857552630417675
random shot 8 : 2.980195250559442
random shot 9 : 2.9117698534023706

shots=[ballistic_sim(normal(vi,Evi),normal(theta,Etheta),normal(yi,Eyi)) for i in range(50)]
print("Range: ",numpy.mean(shots),"+-", numpy.std(shots),"m")
plt.hist(shots,30)
plt.xlabel('Range (m)')
plt.ylabel('Frequency')
plt.title('Range Distribution')

plt.savefig("img/02-range.png",bbox_inches='tight')
plt.close()
#plt.show()

Range:  2.9589400681482583 +- 0.0942704914022101 m

Checkpoint:

• Make sure you can run the above code! (you must use your own data for the abstract)
• Propagate error for the initial velocity for the horizontal shot
• Record values for the angle and height (with measurement error)
• Use Monte Carlo to calculate the expected range for the inclined shot