%matplotlib inline
import numpy as np
import matplotlib.pyplot as pl
#
# Below is the raw data for the position as a function of time for the muffin cup. The "split" function was used to
# separate the data into a list of strings. It then will print out a statement explaining the number of data points.
#
slist=[]
tlist=[]
data = """0.000000000E0 -2.688162330E0
3.336670003E-2 -4.301059729E0
6.673340007E-2 -5.376324661E0
1.001001001E-1 -6.989222059E0
1.334668001E-1 -1.129028179E1
1.668335002E-1 -1.451607658E1
2.002002002E-1 -2.043003371E1
2.335669002E-1 -2.526872591E1
2.669336003E-1 -3.118268303E1
3.003003003E-1 -3.870953756E1
3.336670003E-1 -4.623639208E1
3.670337004E-1 -5.430087907E1
4.004004004E-1 -6.236536606E1
4.337671004E-1 -7.150511799E1
4.671338005E-1 -8.010723744E1
5.005005005E-1 -8.924698937E1
5.338672005E-1 -9.892437376E1
5.672339006E-1 -1.080641257E2
6.006006006E-1 -1.177415101E2
6.339673006E-1 -1.274188945E2
6.673340007E-1 -1.370962788E2
7.007007007E-1 -1.467736632E2
7.340674007E-1 -1.575263126E2
7.674341008E-1 -1.672036969E2
8.008008008E-1 -1.768810813E2
8.341675008E-1 -1.865584657E2
8.675342009E-1 -1.973111150E2
9.009009009E-1 -2.075261319E2
9.342676009E-1 -2.182787812E2
9.676343010E-1 -2.284937981E2
""".splitlines() # split this string on the "newline" character.
print("We have", len(data), "data points.")
We have 30 data points.
#
# This code will take the lists created above and break them into real numbers with reasonable units.
#
tlist = []
ylist = []
for s in data:
t,y = s.split() # break string in two
t=float(t) # convert time to float
y=float(y)/100.0 # convert distanct (in meters) to float
tlist.append(t)
ylist.append(y)
print ("tlist=",tlist)
print ("ylist=",ylist)
tlist= [0.0, 0.03336670003, 0.06673340007, 0.1001001001, 0.1334668001, 0.1668335002, 0.2002002002, 0.2335669002, 0.2669336003, 0.3003003003, 0.3336670003, 0.3670337004, 0.4004004004, 0.4337671004, 0.4671338005, 0.5005005005, 0.5338672005, 0.5672339006, 0.6006006006, 0.6339673006, 0.6673340007, 0.7007007007, 0.7340674007, 0.7674341008, 0.8008008008, 0.8341675008, 0.8675342009, 0.9009009009, 0.9342676009, 0.967634301]
ylist= [-0.0268816233, -0.04301059729, -0.05376324661, -0.06989222059, -0.1129028179, -0.1451607658, -0.2043003371, -0.2526872591, -0.3118268303, -0.3870953756, -0.4623639208, -0.5430087907, -0.6236536606, -0.7150511799, -0.8010723744, -0.8924698937, -0.9892437376, -1.0806412570000001, -1.177415101, -1.274188945, -1.370962788, -1.4677366319999998, -1.575263126, -1.672036969, -1.768810813, -1.865584657, -1.97311115, -2.075261319, -2.182787812, -2.284937981]
#
# This will plot the raw data above. The x-axis will use the tlist and the y-axis will use the ylist.
#
pl.plot(tlist, ylist)
pl.title("raw data")
pl.xlabel("time(s)")
pl.ylabel("y(m)")
pl.grid()
vlist = [] # Velocity list (computed velocities from experimental data)
tvlist = [] # time list (times for corresponding velocities)
for i in range(1,len(tlist)):
dy=ylist[i]-ylist[i-1] # This calculates the change in the ylist.
dt=tlist[i]-tlist[i-1] # This calculates the change in the tlist.
vlist.append(dy/dt) # This addes all the dy/dt values to the vlist
tvlist.append((tlist[i]+tlist[i-1])/2.0)
# This plots the new values of (tv,v)
pl.plot(tvlist,vlist,'g.')
pl.title("Velocity graph")
pl.xlabel("time(s)")
pl.ylabel("$v_y$ (m/s)")
pl.grid()
m=0.0035 # kg of the cupcake
g=9.8 # m/s of gravity
b=0.00375 # drag coefficient
v=0.0 # initial velocity
dt = (tlist[-1]-tlist[0])/(len(tlist)-1) # time per frame in original video
t=0.0
vclist = [v]
tclist = [t]
def deriv(v, t):
return b*v**2/m - g # dv/dt equation
for i in range(len(tlist)): #Creating the Euler's Method graph
dv = deriv(v,t)*dt
v += dv
t += dt
vclist.append(v)
tclist.append(t)
pl.title("Comparison of experimental and drag model")
pl.xlabel("time(s)")
pl.ylabel("velocity (m/s)")
pl.plot(tclist, vclist, 'r-',tvlist,vlist,'g.') #comparing the Euler's method to the raw data
#
# Analytical solution vs Euler Model
#
v=0.0 # initial velocity Euler
v2=0.0 # initial velocity Analytical
dt = (tlist[-1]-tlist[0])/(len(tlist)-1) # time per frame in original video
t=0.0
t2=0.0
vclist = [v]
tclist = [t]
vc2list = [v2] #Analytical vclist
tc2list = [t2] #Analytical tclist
def deriv(v, t):
return b*v**2/m - g # dv/dt equation
def funct(v,t): #I derived this equation by taking the integral of the dv/dt equation. Since everything is a constant
return b*v**2*t/m - g*t # the integral was easy. All I had to do was add a t to the components.
for i in range(len(tlist)): #Creating the Euler's Method graph
dv = deriv(v,t)*dt
vt = funct(v2,t2)
v += dv
t += dt
vclist.append(v)
tclist.append(t)
vc2list.append(v)
tc2list.append(t)
fig, axs = pl.subplots(2) # For this graph setup I used the referenced this website for the code:
axs[0].set_title('Euler Model') # https://matplotlib.org/devdocs/gallery/subplots_axes_and_figures/subplots_demo.html
axs[1].set_title('Analytical Solution') # Since the graphs were so similar, the graphs would be overlapping each other
axs[0].plot(tclist, vclist,'r-') # Because of this, I compared them on top of each other
axs[1].plot(tc2list,vc2list, 'g-')
for ax in axs.flat:
ax.set(xlabel='time(s)', ylabel='velocity(m/s)')
for ax in axs.flat:
ax.label_outer()