P2820 Final Exam W2014

Name:
Student number:

There are 3 questions for a total of 50 marks.  This is a 2-hour, open book exam.

q1 prep (ignore)

Remove["Global`*"]

m=2;
V[t_]=v[t]/.DSolve[{v'[t]==-9.8-0.492v[t]/m,v[0]==5},v[t],t]//First

xtmp=Integrate[V[tt],{tt,0,t}]

X[t_]=xtmp+25

Plot[X[t],{t,0,3},PlotRange->All]

Plot[V[t],{t,0,3},PlotRange->All]

Plot[V'[t],{t,0,3},PlotRange->All]

ParametricPlot[{V[t],V'[t]},{t,0,3}]

vtab=Table[{t,V[t]+0.05 (Random[]-0.5)},{t,0,3,0.1}]

ListPlot[vtab]

SetDirectory["/Users/saika/COURSES/P2820-W2014/P2820FINAL"]

Export["vdata.dat",vtab]

Directory[]

Q1	[17 marks]

Download the file vdata.dat from the course website.  It contains two columns of data: time (s) in the first, velocity (m/s, positive means upwards) in the second.  The data represent a projectile fired straight up.  Acceleration due to gravity is not known.

a) 3 marks

Import and plot the data.  Include axes labels.

Remove["Global`*"]

SetDirectory["<your_directory>"]

FileNames["v*"]

vdata=Import["vdata.dat","Table"];

ListPlot[vdata,AxesLabel->{"time [s]","velocity [m/s]"}]

b) 5 marks

Determine the accelation of the projectile as a function of time using the (discrete form of) the forward difference approximation for the first derivative of a function.  Plot the acceleration as a function of time.

fd[d_,i_]:=((d[[i+1,2]]-d[[i,2]])/(d[[i+1,1]]-d[[i,1]]))

accel=Table[{vdata[[i,1]],fd[vdata,i]},{i,1,Length[vdata]-1}];

ListPlot[accel]

c) 6 marks

Assume that the drag force on the projectile can be approximated with -b v, where v is the velocity and b is the drag coefficient.  Report the drag coefficient b and the gravitational accelaration g to three figures by appropriately plotting the data and determining the slope and intercept through a linear least-squares fit.  Report uncertainties in these quantities. The mass of the projectile is 2kg.

vdata[[1]]

accel[[1]]

m=2;

XY=Table[{vdata[[i,2]],m accel[[i,2]]},{i,1,Length[accel]}];

lp=ListPlot[XY]

f=Fit[XY,{1,x},x]

Show[lp,Plot[f,{x,-20,5}]]

lm=LinearModelFit[XY,{1,x},x]

lm["ParameterTable"]

g=-lm["BestFitParameters"][[1]]/m
lm["ParameterErrors"][[1]]/m

dragcoeff=-lm["BestFitParameters"][[2]]
lm["ParameterErrors"][[2]]

Drag coefficient is 0.489 +/- 0.014  Ns/m
g is 9.69 +/- 0.07  m/s^2

Uncertainties not required.

d) 3 marks

If the initial height of the projectile is 25m, how high off the ground is the projectile at t=3s (the end of the data sampling)?  Use the trapezoid rule to integrate the data.

DNITrap[data_]:=(dt=data[[2,1]]-data[[1,1]];n=Length[data]-1;dt(((data[[1,2]])/2)+((data[[n+1,2]])/2)+Sum[data[[i,2]],{i,2,n}]))

dy=DNITrap[vdata]

yfinal=25+dy

Q2	[24 marks]

A mass of 3.2 kg is attached to a non-ideal spring with potential energy (in J) given by 

U(x) =-1-((95 x)/12)-((73 (x^2))/24)+((17 (x^3))/12)+((13 (x^4))/24), where x is the position of the mass.  Additionally, there is a oscillating external driving force given by,

F(t)= Sin(5 t) as well as a 

drag force of magnitude 0.8 (v^2) (acting in the opposite direction to the motion of the mass), where v is the velocity of the mass.

a) 2 marks

What is the (conservative) force on the particle due to the potential energy?  Write a function that returns this conservative force as a function of the position of the mass.

Remove["Global`*"]

U[x_]:=-1-((95 x)/12)-((73 (x^2))/24)+((17 (x^3))/12)+((13 (x^4))/24)

Fc[x_]=-D[U[x],x]

Plot[U[x],{x,-4,3}]

Plot[Fc[x],{x,-4,3}]

b) 3 marks

What is the net force on the particle?  Write a function that returns the net force as a function of the position of the mass, its velocity and the time.  

Use your function to report the initial force on the particle, given that at t=0 its position is 

xinit = 3 m and its velocity is vinit = -10 m/s.

Fnet[x_,v_,t_]:=Fc[x]+Sin[5 t]- 0.8 v Abs[v]

xinit=3; vinit=-10;tinit=0;

Fnet[xinit,vinit,tinit]

c) 5 marks

Use Mathematica’s built-in ODE solving capability to find the mass’s position as a function of time up to tf = 80s.  Plot the position as a function of time.

tf=80;
m=3.2;

X[t_]=x[t]/.NDSolve[{m x''[t]==Fnet[x[t],x'[t],t],x[0]==xinit,x'[0]==vinit},x[t],{t,0,tf}][[1]]

Plot[X[t],{t,0,80},PlotRange->All]

d) 6 marks

Repeat part d, except now use the 4th-order Runge-Kutta algorithm from class.  Use a timestep of Δt=0.05.

Δt=0.05;
nmax=Round[tf/Δt];
x[1]=xinit;
u[1]=vinit;
t[1]=tinit;

fx[x_,u_,t_]:=u
fu[x_,u_,t_]:=Fnet[x,u,t]/m

Do[t[n+1]=t[n]+Δt;
sx1=fx[x[n],u[n],t[n]];
su1=fu[x[n],u[n],t[n]];
xtmp=x[n]+Δt sx1/2;
utmp=u[n]+Δt su1/2;
sx2=fx[xtmp,utmp,t[n]+Δt/2];
su2=fu[xtmp,utmp,t[n]+Δt/2];
xtmp=x[n]+Δt sx2/2;
utmp=u[n]+Δt su2/2;
sx3=fx[xtmp,utmp,t[n]+Δt/2];
su3=fu[xtmp,utmp,t[n]+Δt/2];
xtmp=x[n]+Δt sx3;
utmp=u[n]+Δt su3;
sx4=fx[xtmp,utmp,t[n]+Δt];
su4=fu[xtmp,utmp,t[n]+Δt];
x[n+1]=x[n]+Δt (sx1+2 sx2+2 sx3+sx4)/6;
u[n+1]=u[n]+Δt(su1+2su2+2su3+su4)/6,
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];

Plot[xN[t],{t,0,tf},PlotRange->All]

e) 4 marks

From your solution from either part c or d, plot the mechanical energy as a function of time.  When does the greatest loss of mechanical energy occur and why is this so?

KE=1/2 m X'[t]^2;
PE=U[X[t]];

Plot[KE+PE,{t,0,tf},PlotRange->All]

Plot[X'[t],{t,0,tf},PlotRange->All]

f) 2 marks

From your solution from either part c or d, determine the period of oscillation near t=70s.  [You can do this by finding the time difference between consecutive peaks.]

Plot[X[t],{t,60,80}]

t2=t/.FindMaximum[X[t],{t,70}][[2]]

t1=t/.FindMaximum[X[t],{t,67.7}][[2]]

period=t2-t1

g) 2 marks

Conisder the potential energy between x=1 and x=2.  What period would you expect if you approximate the oscillating mass as a simple harmonic oscillator?  [Hint: For a harmonic oscillator U = 1/2 k x^2, so k equals U’’(x) evaluated at the minimum.]

Plot[U[x],{x,1,2}]

sol=NSolve[U'[x]==0,x]

xeq=x/.sol[[3]]

U'[xeq]

k=U''[xeq]

2π/(Sqrt[k/m])

Q3	[9 marks]

Consider heat conduction along a straight, insulated bar.  The time evolution of the temperature η along the bar (x is a position along the bar) satisfies a PDE called the heat equation,

		((∂ η)/(∂ t))=α(((∂^2) η)/(∂ (x^2)))

on x ϵ [0, L], t ϵ [0, 15s], where α=0.1(m^2)/s is the thermal diffusivity of the bar, L=3m is the length of the bar, and initial and boundary conditions are given by

	η(x,t=0)=300 Exp[-(((x-1)^2)/0.05)]  (initial temperature profile in degrees celsius)
	
	((∂ η(x=0,t))/(∂ x))=0
	
	((∂ η(x=L,t))/(∂ x))=0.
	
The boundary conditions mean that no heat enters or leaves the bar.

a) 7 marks

Use NDSolve to solve for the temperature of the bar as a function of position and time.

Remove["Global`*"]

c=0.1;
L=3;
T=15;

u[x_,t_]=η[x,t]/.NDSolve[{(∂_t)η[x,t]==c(∂_(x,x))η[x,t],η[x,0]==300 Exp[-(((x-1)^2)/0.05)],Derivative[1,0][η][0,t]==0,Derivative[1,0][η][L,t]==0},η[x,t],{x,0,L},{t,0,T}][[1]]

Manipulate[Plot[u[x,t],{x,0,L},PlotRange->{0,300}],{t,0,T}]

b) 2 marks

Plot the temperature as a function of x at t=15s.  
What is the difference in temperature between the two ends of the bar at t= 15s?

Plot[u[x,T],{x,0,L}]

Δt=u[0,T]-u[L,T]