Mechanics
Mathematica Documentation
Mathematica Source Code
Page 1

During spring training in the middle 1910s, as a publicity stunt, Wilbert Robinson was going to catch a baseball dropped from an airplane. Casey Stengel arranged to surprise Uncle Wilbert by having a grapefruit dropped instead. Take the mass of a grapefruit to be 0.26 kg, and its diameter as 0.12 m. Assume the grapefruit was dropped from a height of 135 m [unrealistically assume it came straight down].

a) Plot its velocity as a function of height assuming no air resistance.
b) Plot its velocity as a function of height assuming an air resistance according to the formula: F(air) = -(1/2)(cross-sectional area)(density of air)(velocity)^2. (The density of air is 1.2 kg/m^3)
c) What is the 'terminal velocity' in m/s expected for grapefruit with air resistance? [Sanity-check your plot in part b) in the light of this 'terminal velocity'.] [The grapefruit missed Robinson's glove, burst on his chest, knocked him down, and made him think for a while that he had been mortally wounded.]

Solution
Part A

Kinematic equation for falling bodies
eq1 = v^2 == v0^2 + 2*a*(x-x0)

Solve the above equation for v
sol = Solve[eq1,v]

If positive x is defined as up, then the velocity of a falling object must be negative. Assign the negative solution.
velocity = -(v0^2 + 2*a*x - 2*a*x0)^(1/2)


Substitute in values into above equation.
velocity1 = velocity/.{a -> -98/10, x0 -> 135, v0 -> 0}
Plot[velocity1,{x,0,135},PlotLabel -> "velocity vs. height"]


Part B

Clear all variables.
ClearAll[eq1,sol,velocity,velocity1,v]

Define constants.
g := 98/10
m := 26/100
dia := 12/100
density := 12/10
Fair := 1/2*crossarea*density*v^2

Find the cross-sectional area.
crossarea = Pi*(dia/2)^2

Newton's Second Law
f = Fair-m*g == -m*a

Change above equation into differential form.
DE = f/.{v -> x'[t],a -> -x''[t]}

Solve the D.E. without initial conditions initial conditions. NOTE: If the initial conditions
are put into the DSolve statement, Mathematica can't solve the DE. Needs["Calculus`DSolve`"]
sol = DSolve[DE,x[t],t]

In order to solve for C[1], the derivative of x[t] must be computed.
dsol = D[sol,t]

Now substitute in the fact that t =0.
dsol0 = dsol/.t -> 0

Assign the solution.
solu = (-910*(3/(26*Pi))^(1/2)* Tanh[(910*(3/(26*Pi))^(1/2)*C[1])/9])/9

Now solve for when solu = 0 in terms of C[1].
solc1 = Solve[solu == 0,C[1]]

Now substitute in the obtained values for C[1] and C[2]. (C[2] = 135)
sol = sol/.{C[1] -> 0, C[2] -> 135}

Assign the solution.
xt = 135 - (3250* Log[Cosh[(21*((3*Pi)/26)^(1/2)*t)/25]])/ (27*Pi)

Velocity is the first derivative of position.
vt = D[xt,t]
ParametricPlot[{xt//N,vt//N},{t,0,8.35},PlotLabel -> "Velocity vs height", PlotRange -> {-20,0}]


Part C

Find the terminal velocity by taking the limit of vt at t approaches infinity. terminalvel = vt/. t -> Infinity//N

At terminal velocity, Fair should equal m*g.
eq2 = Fair == m*g

Solve the above equation for v.
sol = Solve[eq2,v]//N

Thus, the solution checks; the negative answer is the same as the answer to the limit of the velocity as t approaches infinity.

Given Vx(t) = 5.740 cos (1.000 t),
a) Make a graph of average x-acceleration <ax> vs time for a series of time intervals all starting with t0=0.5000 sec. The end of each time interval is t1, which should run between 0.600 and 1.500 sec in steps of 0.1 sec.
b) Why is none of the <ax>> values in part a) equal to the 'true' ax at t=0.5000 sec?
c) make another sequence of <ax> values so that <ax> at t=0.5000 sec will be accurate to four significant figures d) compare your answer in c) to the derivative of V(t), evaluated at t=0.5000 sec.

Solution

Velocity of particle at any time.
v = 574/100*Cos[t]

Equation for average acceleration
accavg = (vf - v0)/(t - t0)

Define constants for this problem.
t0 = 1/2
v0 = v/.t -> 1/2
vf = v

Display contents of accavg
accavg = accavg

Equation for instantaneous acceleration
a = D[v,t]

Time between calculated points on average acceleration graph
timestep = 1/10

Generate 10 points from the average acceleration function to be plotted.
Location of first point
t = t0 + timestep//N
L = {t, accavg}//N

Make a list of ten points starting with the first point.
listofl = Table[{t,accavg},{t,t0+timestep,3/2, timestep}]//N

Generate a graph of the ten points from the average velocity data (but don't display the plot yet).
g := ListPlot[listofl, PlotRange -> {{0,1.6}, {-5.5,-2.5}}]

Generate a graph of the instantaneous acceleration.
h := Plot[a, {t,t0,1.5}, PlotRange -> {{0,1.6}, {-5.5,-2.5}}]

Display both graphs
Show[g,h]

Find minimum time step needed for accavg to be accurate to four significant figures.
Find where accavg and "a" differ by only 0.0005 meters/sec^2.
First, clear the value of the variable t.
ClearAll[t]
sol = FindRoot[accavg - a - 0.0005 == 0, {t,1}] {t -> 0.500199}

Assign the solution
sol = .500199

Calculate the change in time of the solution from t0.
timestep = sol - t0

Find where the tenth point will be placed on the average velocity graph
tmax = 1/2+timestep*10

Generate points as before, but this time with different values of timestep and tmax.
First, assign a value to t
t = t0+timestep 0.500199

listofl = Table[{t,accavg},{t,t0+timestep,tmax, timestep}]//N

Generate a graph of the ten points from the average velocity data (but don't display the plot yet)
g := ListPlot[listofl, PlotRange -> {{.5,.502}, {-2.758, -2.75}}]

Generate a graph of the instantaneous acceleration.
h := Plot[a, {t,t0,tmax}, PlotRange -> {{.5,.502}, {-2.758, -2.75}}]

Display both graphs
Show[g,h]

Values of accavg and "a" at t0, after clearing the value for the variable t
ClearAll[t]
accavg0 = N[accavg/.t -> t0 + timestep,10]
a0 = N[a/.t -> 1/2,10]

Subtract the two to check the number of significant figures
difference = a0 - accavg0
timestep

Thus, accavg0 is accurate to four significant figures when the time step is 0.000199 sec.

Given x(t) = 5.740 sin (1.000 t),
a) Make a graph of average x-velocity <Vx> vs time for a series of time intervals all starting with t0=0.5000 sec. The end of each time interval is t1, which should run between 0.600 and 1.500 sec in steps of 0.1 sec.
b) Why is none of the <Vx>> values in part a) equal to the 'true' Vx at t=0.5000 sec?
c) make another sequence of <Vx> values so that <Vx> at t=0.5000 sec will be accurate to four significant figures d) compare your answer in c) to the derivative of x(t), evaluated at t=0.5000 sec.

Solution

Position of particle at any time
x = 574/100*Sin[1*t]

Equation for average velocity
vavg = (xf - x0)/(t - t0)

Define constants for this problem
t0 = 1/2
x0 = x/.t -> 1/2
xf = x

Display vavg again
vavg

Equation for instantaneous velocity
v = D[x,t]

Time between calculated points on average velocity graph
timestep = 1/10

Generate 10 points from the average velocity function to be plotted
Location of first point
t = t0 + timestep//N
L = {t,vavg}//N

Make a list of ten points starting with the first point
listofl = Table[{t,vavg},{t,t0+timestep,3/2,timestep}]//N

Generate a graph of the ten points from the average velocity data (but don't display the plot yet)
g := ListPlot[listofl, PlotRange -> {{0,1.6}, {0,5.1}}]

Generate a graph of the instantaneous velocity
h := Plot[v, {t,t0,1.5}, PlotRange -> {{0,1.6}, {0,5.1}}]

Display both graphs
Show[{g,h}, PlotLabel -> "Graph of v and vavg vs. time"]

Find minimum time step needed for vavg to be accurate to four significant figures.
Find where vavg and v differ by only 0.0005 meters/sec. (Use the FindRoot command)
First, clear the value of the variable t.
ClearAll[t]
sol = FindRoot[vavg - v - 0.0005 == 0, {t,1}]

Assign the solution
sol = .500363

Calculate the change in time of the solution from t0.
timestep = sol - t0

Find where the tenth point will be placed on the average velocity graph
tmax = 1/2+timestep*10

Generate points as before, but this time with different values of timestep and tmax.
First, assign a value to t
t = t0+timestep 0.500363

listofl = Table[{t,vavg},{t,t0+timestep,tmax,timestep}]//N

Generate a graph of the ten points from the average velocity data (but don't display the plot yet)
g := ListPlot[listofl, PlotRange -> {{.5,.504}, {5.026, 5.038}}]

Generate a graph of the instantaneous velocity
h := Plot[v, {t,t0,tmax}, PlotRange -> {{.5,.504}, {5.026, 5.038}}]

Display both graphs
Show[{g,h}, PlotLabel -> "Graph of v and vavg vs. time"]

Values of vavg and v at t0, after clearing the value for t
ClearAll[t]
vavg0 = N[vavg/.t -> t0 + timestep,10]
v0 = N[v/. t -> 1/2,10]

Subtract the two to check the number of significant figures
difference = v0 - vavg0
timestep

Thus, vavg0 is accurate to 4 significant figures when the time step is 0.000363 sec.

A ball is thrown vertically up . The ball's height as a function of time is given by y(t) = 6+50t-16t2 where t is in seconds and y is in feet.
a) Plot this function to see if it is reasonable.
b) Plot the derivative with respect to time and verify that it is positive while the ball is going up to the highest point and negative thereafter. Is this what you would expect for the vertical component of the velocity?
c) Predict the nature of a graph of the vertical acceleration as a function of time before doing any calculations. Calculate the derivative of the vertical component of the velocity and compare it to your prediction.

Solution

Enter in the height function.
yt = 6 + 50 t - 16 t^2

Plot the function to see if it is a reasonable model.
Plot[yt, {t,0,4}, PlotRange -> {-50,75}, PlotLabel -> yt]

Calculate the derivative of yt.
dyt = D[yt,t]

Plot the derivative to determine the change in the ball's height.
Plot[dyt,{t,0,4}]

Calculate the derivative of the ball's vertical velocity.
ddyt = D[dyt,t]

Here are the components of a supposedly conservative force.
{Fx = xyz^3, Fy = x^2z^3/2, Fz = 3x^2yz^2/2}
Carry out the line integral of F . dl along the following path in the following steps.
a) Start at (x1,y1,z1) and integrate straight to (x2,y1,z1). Note that dl = i dx on this path. This integral is I1.
b) Integrate from (x2,y1,z1) straight to (x2,y2,z1). Call this integral I2.
c) Obtain integrals I3 and I4 by integrating straight from (x2,y2,z1), to (x1,y2,z1), and then straight from there back to (x1,y1,z1).
d) The sum of the four integrals should be zero.[The integral around any closed path should be zero for a conservative force.]
Note: It has no meaning to integrate all 3 components at once along a closed path. Each component of this integral will automatically be zero. (Why?)[Ans: Each component is a scalar potential. Each component integrated around a closed path will vanish. This integration is not a true dot product operation.]

Solution

Enter in the force vector.
Fx = x*y*z^3
Fy = x^2*z^3/2
Fz = 3/2*x^2*y*z^2

Now integrate F over the closed path (x1,y1,z1) to (x2,y1,z1) to (x2,y2,z1) to (x1,y2,z1) and back to (x1,y1,z1).

Integrate from (x1,y1,z1) to (x2,y1,z1).
F1x := Fx/.y -> y1
I1 = Integrate[F1x,{x,x1,x2}]

Integrate from (x2,y1,z1) to (x2,y2,z1).
F2y := Fy/.x -> x2
I2 = Integrate[F2y,{y,y1,y2}]

Integrate from (x2,y2,z1) to (x1,y2,z1).
F3x := Fx/.y -> y2
I3 = Integrate[F3x,{x,x2,x1}]

Integrate from (x1,y2,z1) to (x1,y1,z1).
F4y := Fy/.x -> x1
I4 = Integrate[F4y,{y,y2,y1}]

The sum of the four line integrals should be 0.
I1 + I2 + I3 + I4

The crank of length L connects the rotating wheel of radius R with the mass M which moves in a slot. Let theta = w t, where w is a constant, and determine the maximum tension and maximum compression force in the crank for L=1.0 m and R = 0.5m. [ L must be greater than R. If R>L, the wheel will try to rip the mass out of its slot in a direction perpendicular to the slot.]
a) Develop a formula for the angle phi in terms of R, L, and theta.
b) Take T to be the force exerted by the crank on the mass M. When T is negative, the crank is in tension, and when T is positive, the crank is in compression. Determine maximum and minimum values of T when R = 0.5m and L=1.0m.

Solution
Part A

Relate theta and phi.
f1 = R*Sin[theta] == L*Sin[phi]

Solve for phi.
f2 = Solve[f1,phi]

The corect mathematical solution is:
f2 = ArcSin[R*Sin[theta]/L]

Part B

Get x in terms of R, L, phi, and theta:
f3 = R*Cos[theta] + L*Cos[phi]

Substitute in the value of phi to get position in terms of theta.
f4 = f3/.phi -> ArcSin[(R*Sin[theta])/L]

Write equation in terms of time using theta=omega*t.
f5 = f4/.theta -> omega*t

Take the second derivative of the position function to find the acceleration.
ax = D[f5,{t,2}]

Relate the acceleration of the mass to the force on L.
f7 = M*ax == T*Cos[phi]

Solve the previous equation for force.
f8 = Solve[f7,T]

Assign the solution.
T = -((M*Sec[phi]*(L^2*omega^2*R^2* Cos[omega*t]^2 - L^2*omega^2*R^2*Sin[omega*t]^2 + omega^2*R^4*Sin[omega*t]^4 + L^3*omega^2*R*Cos[omega*t]* ((L^2 - R^2*Sin[omega*t]^2)/L^2)^(1/2) - L*omega^2*R^3*Cos[omega*t]*Sin[omega*t]^2* ((L^2 - R^2*Sin[omega*t]^2)/L^2)^(1/2)))/ (L*(L^2 - R^2*Sin[omega*t]^2)* ((L^2 - R^2*Sin[omega*t]^2)/L^2)^(1/2)))

Write equation f2 in terms of time.
f9 = f2/.theta -> omega*t

Again substitute in the value of phi to get the thesion in terms of time.
f10 = T/.phi -> f9

Simplify the equation.
forcetime = Simplify[f10]
forcetheta = forcetime/.t -> theta/omega

Put in values for a specific case.
ttime = Simplify[forcetime/.{L -> 1, R -> 1/2,
omega -> 2*Pi, M -> 10}]
ttheta = Simplify[forcetheta/.{L -> 1, R -> 1/2,
omega -> 2*Pi, M -> 10}]

Plot time and theta verses force.
Plot[ttime//N, {t,0,2}, PlotLabel -> "Force vs. time"]
Plot[ttheta//N, {theta,0,4*Pi}, PlotLabel -> "Force vs. theta"]

Take the first derivative of the theta function to find the maxes and mins.
dtheta = D[ttheta,theta]

Values of theta at the min
thetamin = FindRoot[dtheta == 0, {theta,.3}]

Value of theta at the max.
thetamax = FindRoot[dtheta == 0, {theta, 2}]

Substitute values of theta back into force equation to get the max and min values.
minforce = N[ttheta/.theta -> 0]
maxforce = N[ttheta/.theta -> 1.92217]

Function to be animated
T1 = forcetheta/.{omega -> 2*Pi, M -> 10}

Function to be plotted in 3D
T2 = T1/.R -> 1/2

Plot3D[T2, {theta,0,2*Pi},{L,3/5,4},PlotLabel -> "Force vs. L and theta", BoxRatios -> {1,1,1}, ViewPoint->{-2.733,-1.587,1.210}]

D1 = Plot3D[T1/.R -> 0, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D2 = Plot3D[T1/.R -> 1/16, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D3 = Plot3D[T1/.R -> 1/8, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D4 = Plot3D[T1/.R -> 3/16, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D5 = Plot3D[T1/.R -> 1/4, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D6 = Plot3D[T1/.R -> 5/16, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D7 = Plot3D[T1/.R -> 3/8, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction ->$DisplayFunction]

D8 = Plot3D[T1/.R -> 7/16, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

D9 = Plot3D[T1/.R -> 1/2, {theta,0,2*Pi},{L,3/5,4}, PlotRange -> {{0,2*Pi},{0,4},{-400,500}}, BoxRatios -> {1,1,1},ViewPoint->{-2.733,-1.587,1.210}, DisplayFunction -> $DisplayFunction]

Now, for an animation of the mass and crank:
Lx = f5
Ly = 0
Rx = R*Cos[omega*t]
Ry = R*Sin[omega*t]

Define the constants.
R := 1/2
L := 1
omega := 2*Pi
M := 10

Here is a simple animation of just the two arms:
Do[Show[Graphics[ Line[{{0,0},{Rx//N,Ry//N},{Lx//N,Ly}}], DisplayFunction -> $DisplayFunction, Prolog -> AbsolutePointSize[4], PlotRange -> {{-1,2},{-1.5,1.5}}, Axes -> True, AspectRatio-> 1]], {t,0,1,1/30}]

The following animation is much more sophisticated.
Do[Show[Graphics[{ Line[{{(R*Cos[omega*t+Pi])//N,(R*Sin[omega*t+Pi]) //N},{Rx//N,Ry//N},{Lx//N,Ly}}], Line[{{(-R*Cos[omega*t+Pi/3])//N,(-R*Sin[omega*t+Pi/3]) //N},{(R*Cos[omega*t+Pi/3])//N,(R*Sin[omega*t+Pi/3]) //N}}], Line[{{(-R*Cos[omega*t+2*Pi/3])//N,(-R*Sin[omega*t+2*Pi/3]) //N},{(R*Cos[omega*t+2*Pi/3])//N,(R*Sin[omega*t+2*Pi/3]) //N}}], Line[{{.3,.1},{1.7,.1},{1.7,-0.1},{.3,-0.1},{.3,.1}}], Circle[{0,0},R], Rectangle[{(Lx-.1)//N,-0.1},{(Lx+.1)//N,0.1}], Rectangle[{1.8,0},{2.0,forcetime/150//N}], Text["Compression", {1.9,1.9}], Line[{{1.8,-300/150},{2.0,-300/150}}], Line[{{1.8,-200/150},{2.0,-200/150}}], Line[{{1.8,-100/150},{2.0,-100/150}}], Line[{{1.8,0},{2.0,0}}], Line[{{1.8,100/150},{2.0,100/150}}], Line[{{1.8,200/150},{2.0,200/150}}], Text["-300 N",{2.1,-300/150},{-1,0}], Text["-200 N",{2.1,-200/150},{-1,0}], Text["-100 N",{2.1,-100/150},{-1,0}], Text["0 N",{2.1,0},{-1,0}], Text["100 N",{2.1,100/150},{-1,0}], Text["200 N",{2.1,200/150},{-1,0}] }], DisplayFunction -> $DisplayFunction, Prolog -> AbsolutePointSize[4], PlotRange -> {{-3.1,3.1},{-3.1,3.1}}, AspectRatio-> 1], {t,0,1,1/30}]

Page 2 of Mathematica Mechanics Problems
Page 3 of Mathematica Mechanics Problems
Page 4 of Mathematica Mechanics Problems