MAPLE text and input lines

```PHYSICS RESOURCE PACKETS PROJECT            File: AIRESIST.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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:
restart;

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.
velocity:=sol[2];

Substitute in values into above equation.
velocity1:=subs({a=-98/10,x0=135,v0=0},velocity);

plot(velocity1,x=0..135,labels=[x,v],title=`velocity vs height`);

Part B

clear all variables.
restart;

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:=subs({v=diff(x(t),t),a=-diff(x(t),t\$2)},f);

Solve the D.E. using initial conditions
sol:=dsolve({DE,x(0)=135,D(x)(0)=0},x(t),numeric);

odeplot(sol,[x(t),diff(x(t),t)],0..8.35,numpoints=150,
title=`Velocity vs height`,labels=[x,v]);

Part C

look at the solution to the DE when t = 1
sol(1);```
```Define velocity as a function of time from the solution to the DE
velocity:=t->rhs(sol(t)[3]);```
```look at the velocity then t = 1
velocity(1);```
```Find the terminal velocity by putting in a large value to t.
terminal_vel:=velocity(10000);```
```At terminal velocity Fair should equal m*g
eq2:=Fair=m*g;```
```Solve the above equation for v
sol:=solve(eq2,v);
evalf(sol[2]);```
```Solution checks!

PHYSICS RESOURCE PACKETS PROJECT            File: AVEACC.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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:
restart;
with(plots):

Velocity of particle at any time.
v:=5.740*cos(1*t);

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

Define constants for this problem.
t0:=1/2;
v0:=subs(t=1/2,v);
vf:=v;

Display contents of accavg
'accavg'=accavg;

Equation for instantaneous acceleration
a:=diff(v,t);

Time between calculated points on average acceleration graph.
time_step:=1/10;

Generate 10 points from the average acceleration function to be plotted.

Location of first point.
t:=t0+time_step:
L:=[t,accavg]:

Make a list of ten points starting with the first point.
for t from (t0+2*time_step) to 1.5 by time_step do
L:=L,[t,accavg];
od:

Display the list.
map(evalf,[L]);

reset the t variable
t:='t':

Generate a graph of the ten points from the average acceleration data.
g:=plot([L],X=t0..1.5,style=point,color=red):

Generate a graph of the instantaneous acceleration.
h:=plot(a,t=t0..1.5):

Display both graphs.
display({g,h},title=`Graph of v and accavg vs time`);

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.
sol:=fsolve(accavg-a=0.0005,t=t0..1.5);

Calculate the change in time of the solution from t0.
time_step:=sol-t0;

Find where the tenth point will be placed on the average acceleration graph
tmax:=1/2+time_step*10;

Generate points as before, but this time with different values of time_step and
tmax.
t:=t0+time_step:
L:=[t,accavg]:
for t from (t0+2*time_step) to tmax by time_step do
L:=L,[t,accavg];
od:
map(evalf,[L]);
t:='t':

g:=plot([L],Time=t0..tmax,style=point,color=red):
h:=plot(a,t=t0..tmax):

display({g,h},title=`Graph of v and accavg vs time`);

Values of accavg and "a" at t0

L[1][2] (second item of the first item in list L; the first value of accavg)
accavg0:=evalf(L[1][2]);
a0:=evalf(subs(t=1/2,a));

Subtract the two to check the number of significant figures
difference=a0-accavg0;
TimeStep=time_step;

accavg0 is thus accurate to 4 significant figures when the time step is
0.0001985325 sec.

PHYSICS RESOURCE PACKETS PROJECT            File: AVEVEL.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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:
restart;
with(plots):
Position of particle at any time.
x:=5.740*sin(1*t);

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

Define constants for this problem.
t0:=1/2;
x0:=subs(t=1/2,x);
xf:=x;

Display contents of vavg
'vavg'=vavg;

Equation for instantaneous velocity
v:=diff(x,t);

Time between calculated points on average velocity graph.
time_step:=1/10;

Generate 10 points from the average velocity function to be plotted.

Location of first point.
t:=t0+time_step:
L:=[t,vavg]:

Make a list of ten points starting with the first point.
for t from (t0+2*time_step) to 1.5 by time_step do
L:=L,[t,vavg];
od:

Display the list.
map(evalf,[L]);

reset the t variable
t:='t':

Generate a graph of the ten points from the average velocity data.
g:=plot([L],X=t0..1.5,style=point,color=red):

Generate a graph of the instantaneous velocity.
h:=plot(v,t=t0..1.5):

Display both graphs.
display({g,h},title=`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.
sol:=fsolve(vavg-v=0.0005,t=t0..1.5);

Calculate the change in time of the solution from t0.
time_step:=sol-t0;

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

Generate points as before, but this time with different values of time_step and
tmax.
t:=t0+time_step:
L:=[t,vavg]:
for t from (t0+2*time_step) to tmax by time_step do
L:=L,[t,vavg];
od:
map(evalf,[L]);
t:='t':

g:=plot([L],Time=t0..tmax,style=point,color=red):
h:=plot(v,t=t0..tmax):

display({g,h},title=`Graph of v and vavg vs time`);

Values of vavg and v at t0

L[1][2] (second item of the first item in list L; the first value of vavg)
vavg0:=evalf(L[1][2]);
v0:=evalf(subs(t=1/2,v));

Subtract the two to check the number of significant figures
difference=v0-vavg0;
TimeStep=time_step;

vavg0 is thus accurate to 4 significant figures when the time step is
0.000363224 sec.

PHYSICS RESOURCE PACKETS PROJECT            File: BALLUP.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A ball is thrown vertically up.  The ball's height as a function of time is
given by  y(t) = 6+50t-16t^2  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:
restart;
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);

Calculate the derivative of yt.
dyt:= diff(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:=diff(dyt,t);

PHYSICS RESOURCE PACKETS PROJECT            File: CONSHORT.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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.]

Solution:
restart;
with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace

Enter in the force vector.
F:=vector([x*y*z^3,x^2*z^3/2,3*x^2*y*z^2/2]);

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

Force on particle in +x direction from A to B
F1:=subs({y=y1,z=z1},F[1]);

Work done on particle from A to B
I1:=int(F1,x=x1..x2);

Force on particle in +y direction from B to C
F2:=subs({x=x2,z=z1},F[2]);

Work done on particle form B to C
I2:=int(F2,y=y1..y2);

Force on particle in x direction from C to D
F3:=subs({y=y2,z=z1},F[1]);

Work done on particle from C to D
I3:=int(F3,x=x2..x1);

Force on particle in y direction from D to A
F4:=subs({x=x1,z=z1},F[2]);

Work done on particle from D to A
I4:=int(F4,y=y2..y1);

Work around the closed path is 0.
Isum:=I1+I2+I3+I4;

PHYSICS RESOURCE PACKETS PROJECT            File: CRANK.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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:
restart;
with(plots):

Part A

Relate theta and phi.
f1:=R*sin(theta)=L*sin(phi);

Solve for phi.
f2:=phi=solve(f1,phi);

Part B

Get x in terms of R, L, phi, theta:
f3:=x=R*cos(theta)+L*cos(phi);

Substitute in the value of phi to get the position in terms of theta.
f4:=subs(f2,f3);

Simplify the result.
f4:=simplify(f4);

Write equation in terms of time using theta=omega*t.
f5:=subs(theta=omega*t,f4);

Take the second derivative of the position function to find the acceleration.
a[x]:=diff(rhs(f5),t\$2);

Relate the acceleration of the mass to the force on L.
f7:=M*a[x]=T*cos(phi);

Solve the previous equation for force T.
f8:=T=solve(f7,T);

Write equation f2 in terms of time.
f9:=subs(theta=omega*t,f2);

Again substitute in the value of phi to get the force T in terms of time.
f10:=subs(f9,f8);

Simplify the equation.
T_time:=simplify(f10);

T_theta:=subs({t=theta/omega},T_time);

Put in values for a specific case.
sub:=subs({L=1,R=1/2,omega=2*Pi,M=10},[T_time,T_theta]);

Equations for T.
Tt:=simplify(rhs(sub[1]));
Ttheta:=simplify(rhs(sub[2]));

plot(Tt,t=0..2,title=`T vs time`);

plot(Ttheta,theta=0..4*Pi,title=`T vs theta`);

Take first derivative of T to find the maxes and mins.
dT:=simplify(diff(Ttheta,theta));

Value of theta at the min.
thetamin:=fsolve(dT=0,theta=-1..1);

Value of theta at the max.
thetamax:=fsolve(dT=0,theta=1.8..2.2);

Substitute values of theta back into the T force equation to get the max and min
values.
MinT:=evalf(subs(theta=thetamin, Ttheta));
MaxT:=evalf(subs(theta=thetamax, Ttheta));

Function to be animated.
T1:=rhs(subs({omega=2*Pi,M=10},T_theta));

Function to be plotted in 3D.
T2:=subs(R=1/2,T1);

plot3d(T2,theta=0..2*Pi,L=3/5..4,title=`T vs L and theta`);

Animate the 3D plot for different values of R.
animate3d(T1,theta=0..2*Pi,L=3/5..4,R=0..1/2,title=`T vs L, theta, and R`);

Animation of the crank

Equation for the position of the mass.
f5;

Substitute in values for problem.

x and y coordinates of wheel.

Generate frames for the animation.
for t from 0 to 29/30 by 1/30 do

g.(30*t):=plot({[[-Rx,-Ry],[Rx,Ry],[positionM,0]],[[-Rx2,-Ry2],[Rx2,Ry2]],
[[-Rx3,-Ry3],[Rx3,Ry3]],[lengthR*cos(theta),lengthR*sin(theta),theta=0..2*Pi]}
,color=black,thickness=2);
h.(30*t):=polygonplot([[positionM-0.1,0.07],[positionM+0.1,0.07],[positionM+0.1,
-0.07],[positionM-0.1,-0.07],[positionM-0.1,0.07]],color=blue);
i.(30*t):=polygonplot({[[positionM-0.1,0.025],[-0.05,0.025],[-0.05,-0.025],
[positionM-0.1,-0.025]],[[positionM+0.1,0.025],[lengthL+lengthR+0.2,0.025],
[lengthL+lengthR+0.2,-0.025],[positionM+0.1,-0.025]]},color=grey);

crank.(30*t):=display({g.(30*t),h.(30*t),i.(30*t)});
od:
t:='t';

Display the animation.
display([crank.(0..29)],insequence=true,scaling=constrained,axes=frame);

PHYSICS RESOURCE PACKETS PROJECT            File: HAMMER.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

The force imparted by a hammer to a steel bar is given by 50 N (1 - (t/to)^3),
where t runs from 0 to to. [The force is assumed to have suddenly reached 50 N
as the hammer initially struck the bar.] Take to=20 ms and determine the impulse
graphically which is associated with the collision, and also use maple or
mathematica to find the impulse.

Solution:
restart;

Equation for force on steel ball.
f:=50*(1-(t/t0)^3);

Substitute in values into above equation.
F:=subs(t0=20/1000,f);

Find where force equals 0
sol:=solve(F,t);

Assign to tf.
tf:=sol[1];

Plot the equation.
plot(F,t=0..tf,title=`Force on stell ball vs time.`);

Show integral in symbolic form.
Impulse:=Int(F,t=0..tf);

Evaluate integral.
Impulse:=int(F,t=0..tf);

PHYSICS RESOURCE PACKETS PROJECT            File: JIM&PHIL.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

Jim throws a baseball at an angle theta1 with a speed v1. Phil, who is standing
in the vertical plane of motion, also throws a baseball in the same forward
direction at an angle theta2 with a speed v2. Phil throws his ball at the exact
moment Jim's ball is directly overhead. At that moment, Phil notices that Jim's
ball is neither rising nor falling, being at its maximum height.
a) With what speed (as a function of the initial launch angles and v1) must
Phil's ball be thrown to collide with Jim's?
b) Determine where the collision takes place.
c) After the first ball is thrown, how long does it take the two balls to
collide? Assume Jim and Phil are at the same height on a level field.

Solution:
restart;
DETERMINATION OF SPEED FOR A COLLISION (PART A)
Enter in the initial equations given
x1 := v0x1*t1;
x2 := x02 + v0x2*t2;
t1 := t2 + Dt;
where Dt is Delta*t
Start eliminating variables; realize what x02 equals
x02 := v0x1*Dt;

In order for a collision to occur, x1  must equal x2
equat1 := x1=x2;

An expansion will "clean up" the equation
equat2 := expand(equat1);

Now solve for v0x1
equat3 := v0x1=solve(equat2, v0x1);

Now generalize the problem; use v*cos(theta) instead of vx
equat4 := v01*cos(theta1) = v02*cos(theta2);

Solving for v02 gives:
equat5 := v02 = solve(equat4,v02);

DETERMINATION OF LOCATION OF COLLISION (PART B)
Typing in the word 'restart' will wipe out all of the previously defined
variables.
restart;
In order for a collision to occur, both y1 = y2 and x1 = x2 must be true.
equat6 := y1 = y2;

Write the standard expressions for y1 and y2
equat7 := y01 + v0y1*t1 - (g/2)*(t1)^2 = y02 + v0y2*t2 - (g/2)*(t2)^2;

Since both balls were thrown from the same height, y01 and y02 can be canceled
out by setting them equal to 0.
y01 := 0;
y02 := 0;

Define what t1 is again (Maple doesn't remember because of the 'restart'
command)
t1 := t2 + Dt;

Dt can be determined by thinking about the relationship between x02 and the
velocity fo the ball at that point.
Dt := v0y1/g;

Now substitute these terms in to equation 7, the one that sets the 2
y-components equal
equat8 := simplify(equat7);

Solving for the term v0y2*t2 yields:
equat9 :=v0y2* t2 = solve(equat8,v0y2*t2);

Thus, t2 can be determined
equat10 := equat9/v0y2;

Now the location of the collision can be determined (Find when x1 = x2 and
y1 = y2)
equat11 := x2 = x02 + v0x2*t2;

Redefine what x02 and v0x2 are (Look at the solution to part A)
x02 := v0x1*Dt;
v0x2 := v0x1;

Substituting in these values, along with t2, gives:
equat12 := subs(t2 = rhs(equat10),equat11);

'simplify' will help make the equation look better
equat13 := simplify(equat12);

Now for the y-components:
equat14 := y2 = y0 + rhs(equat7);

Plugging in t2 gives:
equat15 := subs(t2 = rhs(equat10),equat14);

where y0 must be specified

TIME FOR PROJECTILES TO COLLIDE (PART C)
Instead of a 'restart' which wipes out all of the variables, the tick marks
allow for a single variable to be redefined to itself.
t1 := 't1';
Dt := 'Dt';
equat16 := t1 = t2 + Dt;

Defining what t2 and Dt are gives:
t2 := rhs(equat10);
Dt := v0y1/g;
equat14 := t1 = t2 + Dt;

'factor' will help "clean up" the final answer
equat15 := factor(equat14);

PHYSICS RESOURCE PACKETS PROJECT            File: LAUNCH.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A ball is launched with an initial velocity of 130 m/s at ground level. It
strikes a wall 150 m away, at a height of 10m.
a) Find the parametric equations y(t) and x(t) in terms of  t and launch angle
theta.
b) Find the launch angle and time of flight. [Note: There are two solutions.]
c) Determine the path y(x) from the parametric equations.
d) Do a 'parametric plot' of x(t) and y(t). Plot y(x) and compare to the
parametric plot.

Solution:
restart;

Kinematic projectile motion equations.
f1:=x=v0*cos(theta)*t;
f2:=y=v0*sin(theta)*t-49/10*t^2;

Substitute in specific values.
xcomp:=subs({v0=130,x=150},f1);
ycomp:=subs({v0=130,y=10},f2);

Solve xcomp and ycomp for t and theta.
sol:=solve({xcomp,ycomp},{t,theta});

Find the all of the solutions to the equations.
sol2:=evalf(allvalues(sol));

Through out values that don't make sense because t<0.
angle1:=theta=subs(sol2[1],theta);
angle2:=theta=subs(sol2[3],theta);

Substitute values of theta into the kinematic equations.
x1:=rhs(subs(angle1,xcomp));
y1:=rhs(subs(angle1,ycomp));
x2:=rhs(subs(angle2,xcomp));
y2:=rhs(subs(angle2,ycomp));

Plot {x1,y1} and {x2,y2}
plot({[x1,y1,t=0..30],[x2,y2,t=0..30]},x=0..150,y=0..900,title=`Parametric
plot of projectile motion`);

Solve first kinematic equation for time.
Time:=solve(f1,t);

Substitute into second kinematic equation.
path:=subs(t=Time,f2);

Substitute angle1 and angle2 into path.
path1:=evalf(subs({v0=130,angle1},path));
path2:=evalf(subs({v0=130,angle2},path));

plot({rhs(path1),rhs(path2)},x=0..150,title=`Paths of projectile as a function
of x`);

Plots are the same.

PHYSICS RESOURCE PACKETS PROJECT            File: MASSINC.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
moloney@nextwork.rose-hulman.edu            Software: Maple V Release 4

Mass M is released from rest as shown in the figure and moves down the plane.
a) Assuming no friction anywhere in the system, how far down the incline will M
move before coming momentarily to rest? [This can be done just using
conservation of energy.]
b) Write down Newton's second law for the mass M on the incline.
c) Solve this equation using 'dsolve' for the general equation of displacement
vs. time.
[In maple, use the commands assume(k>0) and assume(M>0) before executing
the dsolve command]
For parts d) and e), let M=2.5 kg, k=175 N/m, and g=9.8 m/s^2.
d) Plot the displacement vs time for theta = 30 degrees.
e) Plot the displacement vs. time for various angles values of theta from 0 to
90 degrees in increments of 15 degrees. [Hint: A loop would be helpful.]

Solution:
restart;
with(plots):

Part A

Conservation of Energy
energy:=m*g*x*sin(theta)=1/2*k*x^2;

Solve the above equation for x to find the maximum displacement of the mass.
sol:=solve(energy,x);

We're interested in only the second solution.
maxdist:=sol[2];

Maximum displacement of mass at 30deg.
maxdist_at_30_deg:=evalf(subs(m=5/2,k=175,g=98/10,theta=30*Pi/180,maxdist));

Part B

Differential equation describing the motion of the mass.
f:=m*g*sin(theta)-k*x(t)=m*diff(x(t),t\$2);

Part C

Tell Maple that that the mass and the spring constant are positive.
assume(m>0);
assume(k>0);

Solve the differential equation.
sol:=dsolve({f,x(0)=0,D(x)(0)=0},x(t));

Part C

Substitute in values given in problem.
xt:=subs(m=5/2,k=175,g=98/10,rhs(sol));
xt:=collect(xt,sin(theta));

plot(subs(theta=30*Pi/180,xt),t=0..4,title=`Path of mass vs time for angle=30
deg`);

Part D

Generate graphs for the motion of the mass for theta = 0,15,30,45,60,75, and 90
deg.
for th from 0 to 6 do
h.th:=plot(subs(theta=th*15*Pi/180,xt),t=0..4);
od:

display([h.(0..6)],title=`Path of mass vs time for angle=0 deg to 90 deg`);

Find maximum displacement of mass when theta = 30deg.

Maximum displacement occurs when the cosine term of the motion = 1.
maxdisp:=evalf(subs(t*sqrt(70)=Pi,motion30));

The above answer agrees with the one generated using conservation of energy.

plot3d(xt,t=0..2,theta=0..Pi/2,title=`Path of mass vs time for angle=0 deg
to 90 deg`);

PHYSICS RESOURCE PACKETS PROJECT            File: MASSINFR.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
moloney@nextwork.rose-hulman.edu            Software: Maple V Release 4

This is the MASSINC problem with static and kinetic friction coefficients given,
and with the same static and kinetic coefficients.

Let M=2.5 kg, k=175 N/m, theta = 60 degrees, and g=9.8 m/s^2.
a) Write down Newton's second law for the mass M on the incline.
b) Solve this equation using a numeric 'dsolve' for  displacement vs. time.
c) Plot displacement vs time for the mass.

Solution:
restart;

Needed to plot a DE evaluated numerically.
with(plots):
Part A

DE for motion of the mass if there is friction on the plane.
F:=m*g*sin(theta)-k*x(t)-mu*g*cos(theta)*Heaviside(diff(x(t),t))+mu*g
*cos(theta)*Heaviside(-diff(x(t),t))=m*diff(x(t),t\$2);

Substitute in values from problem.
Ft:=simplify(subs(m=5/2,k=175,g=98/10,theta=60*Pi/180,mu=1/4,F));

Part B

Solve the DE numerically
sol2:=dsolve({Ft,x(0)=0,D(x)(0)=0},x(t),numeric);

Part C

odeplot(sol2,[t,x(t)],0..3.4,numpoints=150,title=`position of block vs time
with friction on plane`);

PHYSICS RESOURCE PACKETS PROJECT            File: MSPRING.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A mass M hangs from the ceiling by a  massless spring of constant k and
unstretched length Lo. The length of the spring as a function of time is h(t),
and the extension of the spring is h(t) - Lo.
a)  Show that when the mass is hanging at rest, the spring length is  Lo +Mg/k.
The mass is thus at a distance from the ceiling of h(0) =  Lo + Mg/k.  This
is the 'equilibrium' position of the mass.
When the mass is moving, h(t) = Lo + Mg/k +x(t),  where x(t) is the distance
from the equilibrium position.
b) Let M = 1.55 kg, k = 350 N/m, and Lo = 0.37 m, and let the mass be released
at t=0 when the spring is unstretched. (Note that x(t) will be negative when
this happens.)  Solve Newton's second law applied to the mass M (using
dsolve, or a Runge-Kutta integration) and plot the motion of the mass.
c) Compare the frequency of the motion to the theoretical value of
(k/M)^(1/2) / (2 pi).  (You should also be able to notice that the amplitude
of the motion is Mg/k.)

Solution:
restart;

Part A: Equilibrium Position

Downward pull due to gravity equals the upward pull due to the spring.
eq1:= m*g=k*x;

Solve eq1 for x to find out how much the mass stretches the spring.
stretch:=solve(eq1,x);

Add this to the original length of the spring to find the spring's total
length.
h0:=L[0]+stretch;

Value of h when mass is moving.  x(t) is the distance from the equilibrium
point.
ht:=h0+xt;

Part B: Motion of mass

Newton's second law.
f:=m*a=mg-k*(ht-L[0]);

Change to differential form.
DE:=m*diff(x(t),t\$2)=m*g-k*(m*g/k+x(t));

Assign values for the constants given in the problem.
m:=155/100; k:=350; g:=98/10; L[0]:=37/100;

Now solve the D.E. using the initial conditions.
x(0)= initial position; D(x)(0)= initial velocity
sol:=dsolve({DE,x(0)=-m*g/k,D(x)(0)=0},x(t));

Store the solution in xt.
xt:=rhs(sol);

plot(xt,t=0..2,title=`Motion of mass vs time`);

Part C: Comparison to theoretical value.
theoretical_freq:=combine((k/m)^(1/2)/(2*Pi));

Display in fractional form.
evalf(theoretical_freq);

sol;

Looking at the function for the position of the mass, we see that is frequency
is determined by the cosine term.

The angular velocity is thus:
omega:=10/31*sqrt(2170);

This make the frequency of the function:
freq:=omega/(2*Pi);

Or in decimal form:
evalf(freq);

The frequency of motion from the derived D.E. is the same as the theoretical
value.

Find the value of m*g/k.
m*g/k;

The value of m*g/k is the same as the coefficient in front of the cosine term of
x(t), the amplitude of the oscillation.

PHYSICS RESOURCE PACKETS PROJECT            File: RCSPRING.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

Drive a mass-spring system at various frequencies and observe the frequency
response.  Use a mass of 1 kg, a spring constant of 400 N/m, a damping force of
10 N*s/m, and sin(2*Pi*f*t) for the forcing function.

a) Write the differential equation describing the motion of the mass.  Let x=0
be the equilibrium position of the mass with x being positive in the upward
direction.
b) Plot the steady state frequency response for f = 0 Hz to 10 Hz

Solution:
restart;
with(plots):
assume(f>0);

Part A

Mass on spring (kg)
m:=1;

Gravitational constant
g:=981/100;

Spring constant (N/m)
k:=400;

Force spring exerts on mass (N)
Fs:=-k*x(t);

Damping (N*s/m)
b:=10;

Forcing function on mass
Ff:=sin(2*Pi*f*t);

Newton's second law
DE:=Fs+Ff-b*diff(x(t),t)=m*diff(x(t),t\$2);

Part B

Solve the DE for x(t)
sol:=rhs(dsolve({DE,x(0)=0,D(x)(0)=0},x(t)));

Simplify the result.
xt:=combine(sol,trig);

Find the range of xt for the steady state response.  The format of the maple
output is "Range := min x value..max x value".
Range:=limit(xt,t=infinity);

Pick off the maximum value from Range to find the magnitude of the steady state
frequency response.
mag:=convert(mag,list)[2];

plot(mag1,f=0..10,title=`Frequency Response of the System (Amplitude vs f)`);

PHYSICS RESOURCE PACKETS PROJECT            File: RODGRAV.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A very narrow rod has length L=0.90 m, and mass M=1.75 kg.  Determine the total
gravitational force that this rod exerts on a tiny sphere of mass m=0.01 kg,
located at a distance h=0.55 m from the center of the rod along a perpendicular
bisector.  Do this by finding the force exerted on the sphere by a small segment
of the rod Dx=0.01 m, and then adding all forces from small segments to get the
total force.

Solution:

restart;

Define constants:
L:=0.9;
M:=1.75;
m:=0.01;
h:=0.55;
G:=6.67*10^(-11);

Mass of one segment of the bar.
Ms:=M/9;

Distance from bar segment to sphere mass.
r:=sqrt(x^2+0.55^2);

Force between segment of bar and sphere mass.
F:=G*m*Ms/(r^2);

Force between segment of bar and sphere mass in the y direction.
Fy:=F*h/r;

Add up the force in the y direction due to the nine segments of the rod.

Set the total to zero.
Ftot:=0;

Add the forces from each segment to the total.
for x from -0.4 to 0.4 by 0.1 do
Ftot:=Ftot+Fy;
od;

The total force on the sphere mass is:
Ftot;

PHYSICS RESOURCE PACKETS PROJECT            File: SKBOARD.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

Determine the angle theta that will in principle enable a skate-boarder who
starts at the top (height H) with zero velocity to achieve a maximum jump
distance. [Treat the height H as a constant, along with gravity.] The jump
distance is measured horizontally from the end of the ramp, which has radius R
and extends through and angle   from the vertical.  The height h of the end of
the ramp above the ground is R(1 - cos (theta)).

a) Obtain an expression for the horizontal jump distance L to reach the ground.
[Hint: find the time to reach the ground.]
For parts b) and c)assume  H = 5.1 m.
b) Determine the ramp angle theta for maximum jump distance when R = 1.1 m.
Going Further: (optional)
c) What happens to the maximum jump distance as R decreases?
[Ans: Jump max dist increases as R->0, because it doesn't lose any energy
reaching the launch angle]

Solution:
restart;
with(plots):

Define constants.
g:=98/10;

Height at the end of the ramp
h:=R*(1-cos(theta));

Conservation of energy:
eq1:=m*v[0]^2/2=m*g*(H-h);

Solve the above equation for v0
sol:=solve(eq1,v[0]);

Pick the positive solution.
v[0]:=sol[1];

x and y components of velocity:
v[x]:=v[0]*cos(theta);
v[y]:=v[0]*sin(theta);

Equation for y position of skate-boarder.
eq2:=h+v[y]*t-g*t^2/2;

Solve above equation for time.
sol2:=solve(eq2=0,t);

Pick result with positive value for time.
t:=simplify(sol2[1]);

Part A

Distance traveled by skate-boarder
distance:=v[x]*t;

Substitute in values for H and R into the distance equation.
distance1:=simplify(subs({H=51/10,R=11/10},distance));

plot(distance1,theta=0..Pi/2,title=`Distance traveled vs theta`);

Derivative of the distance equation.
d_dist:=diff(distance1,theta);

Part B

Find theta where the distance traveled is a maximum.
thetamax:=fsolve(d_dist=0,theta=0..Pi/2);

Convert to degrees
thatamax_deg:=evalf(thetamax*180/Pi);

Now substitute only H into the distance equation so that theta and R may be
varied.
distance2:=simplify(subs(H=51/10,distance));

Make a 3d plot of distance2
plot3d(distance2,theta=0..Pi/2,R=0..10,grid=[40,40],title=`Distance vs theta &
R`);

Animate distance2 as R decreases to find out what happens to the jump distance.

The ranges on the animate command must be in increasing order, so to get it to
animate the function as R decreases make a new function by substituting -R for
R.
distance3:=subs(R=-R,distance2);

The plot range for R is now -10 to 0, but this represents 10 to 0 for R.
animate(distance3,theta=0..Pi/2,R=-10..0,labels=[theta,Dist],title=`Distance
vs theta as R decreases`);

The maximum jump distance increases as R decreases.

PHYSICS RESOURCE PACKETS PROJECT            File: TEACUP.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

Certain amusement park rides consist of two circular motions superimposed as
shown in the sketch. There is an arm of length R1 and a second arm of length R2.
The first arm rotates at angular speed w1 and the second arm rotates with
respect to the first arm at an angular speed of w2.

a) Write an expression for the position of point P (at the end of the second
arm) as a function of time.
b) Using values R1 = 10.0 m, R2 = 5.00 m, periods of T1 = 15.0 s, and
T2 = 8.0 s,  and create a parametric (x,y) plot of the motion of point P.
c) Write an expression for the velocity vector of point P as a function of time.
Graph the magnitude of the velocity as a function of the time. Make the time
interval for the plot long enough to view the entire range of possible
velocities. What is the physical configuration of the ride at moments of
greatest velocity?
d) Obtain the general expression for the acceleration vector as a function of
time. Graph the magnitude of the acceleration as a function of time showing
the complete range of values. What is the physical configuration of the ride
at moments of greatest acceleration?
Going further: e) Change the values of the radii and rotation rates to give
customers a ride with a maximum acceleration of 3 g.
Still further: f) Rework parts a) - d) with a negative w2, that is w2 in the
opposite direction of w1.
[Enhancements in solution: 3D plot of velocity magnitude, and acceleration vs x,
y position.]

Solution:
restart;
with(plots):
Position of point Q in general case.
QxGen:=R1*cos(omega[1]*t);
QyGen:=R1*sin(omega[1]*t);

Position of point P with respect to Q in general case.
PrQxGen:=R2*cos((omega[1]+omega[2])*t);
PrQyGen:=R2*sin((omega[1]+omega[2])*t);

Part A

Position of point P

General case
PxGen:=QxGen+PrQxGen;
PyGen:=QyGen+PrQyGen;

Part B

Substitute in values for specific case
values:={R1=10,R2=5,omega[1]=2*Pi/15,omega[2]=2*Pi/8};
Qx:=subs(values,QxGen);
Qy:=subs(values,QyGen);
Px:=subs(values,PxGen);
Py:=subs(values,PyGen);

plot([Px,Py,t=0..120],title=`Parametric plot of motion from T=0-120s`);

Create frames for an animation of the problem. (This may take some time)

(Changing the range to 0-120 will draw the complete cycle of motion, but will
take up more memory and time.)
for t from 0 to 60 by 1/3 do

Value of constant in g.(t*3) is equal to 1/(the time step).  In this case it is
1/(1/3)=3
g.(t*3):=plot([[0,0],[Qx,Qy],[Px,Py]]):
od:

Display the animation.

(The range in g.(0..180) is equal to the total number of frames.  In this case
it is 60/(1/3)=180).
display([g.(0..180)], insequence=true);

Reset t variable
t:='t';

Part C

Find the velocity and speed of point P
vPx:=diff(Px,t);
vPy:=diff(Py,t);
speed:=sqrt(vPx^2+vPy^2);
speed:=simplify(speed);

plot(speed,t=0..50,title=`Speed of P vs time`);

Find where the slope of the curve equals 0.
dspeed:=simplify(diff(speed,t));

Fsolve finds a solution within the interval specified.
Time:=fsolve(dspeed=0,t=15..19);

Location of P at maximum speed.
locationP:=simplify(subs(t=Time,[Px,Py]));
locationQ:=simplify(subs(t=Time,[Qx,Qy]));

Convet to polar coordinates.
locP:=polar(locationP[1]+locationP[2]*I);
locQ:=polar(locationQ[1]+locationQ[2]*I);

Fastest speed is when P is furthest from the center.
plot speed on z axis and position on x and y axes.
spacecurve([Px,Py,speed],t=0..120,numpoints=500,title=`Speed of P vs
position`);

Part D

Find the acceleration of point P
aPx:=diff(vPx,t);
aPy:=diff(vPy,t);
MagAcc:=sqrt(aPx^2+aPy^2);
MagAcc:=simplify(MagAcc);

plot(MagAcc,t=0..50,title=`|a| vs time`);

dMagAcc:=simplify(diff(MagAcc,t));

Fsolve finds a solution within the interval specified.
Time:=fsolve(dMagAcc=0,t=15..19);

Location of P at maximum acceleration.
locationP:=simplify(subs(t=Time,[Px,Py]));
locationQ:=simplify(subs(t=Time,[Qx,Qy]));

Convet to polar coordinates.
locP:=polar(locationP[1]+locationP[2]*I);
locQ:=polar(locationQ[1]+locationQ[2]*I);

Greatest acceleration is when P is furthest from the center.
plot |a| on z axis and position on x and y axes.
spacecurve([Px,Py,MagAcc],t=0..120,numpoints=500,title=`|a| of P vs
position`);

Part E

Change equations back to general case.

"t\$2" takes second derivative
aPx:=diff(PxGen,t\$2);
aPy:=diff(PyGen,t\$2);

General expression for |a|
MagAcc:=sqrt(aPx^2+aPy^2);

We want the maximum |a| to be 3g and there are four variables that effect
acceleration.  Pick three of these variables, and have Maple solve for the
fourth one.

Substituting in values for w1, w2 and R2:
MagAcc1:=subs({omega[1]=2*Pi/8,omega[2]=2*Pi/6,R2=6},MagAcc);
MagAcc1:=simplify(MagAcc1);

Tell Maple we are looking for a positive value of R1.
assume(R1>0);

We know from the previous steps that acceleration is a maximum when t=0, so
substitute t=0 into MagAcc1.
MagAcc2:=simplify(subs(t=0,MagAcc1));

Now find value of R1 that will result in a maximum acceleration of 3g.
sol:=R1=solve(MagAcc2=3*98/10,R1);
soln:=evalf(sol);

Substituting in the value for R1, the equation for |a| as a function of time
becomes:
MagAcc3:=simplify(subs(sol,MagAcc1));

The g force is then:
Gforce:=MagAcc3/(98/10);

Plot the g force
plot(Gforce,t=0..50,title=`g force of ride`);

Repeat steps (a)-(d) for w2 going in the opposite direction.
restart;
with(plots):

Position of point Q in general case.
QxGen:=R1*cos(omega[1]*t);
QyGen:=R1*sin(omega[1]*t);

Position of point P with respect to Q in general case.  (This time subtract
omega2)
PrQxGen:=R2*cos((omega[1]-omega[2])*t);
PrQyGen:=R2*sin((omega[1]-omega[2])*t);

Part A

Position of point P

General case
PxGen:=QxGen+PrQxGen;
PyGen:=QyGen+PrQyGen;

Substitute in values for specific case
values:={R1=10,R2=5,omega[1]=2*Pi/15,omega[2]=2*Pi/8};
Qx:=subs(values,QxGen);
Qy:=subs(values,QyGen);
Px:=subs(values,PxGen);
Py:=subs(values,PyGen);

plot([Px,Py,t=0..120],title=`Parametric plot of motion form T=0-120s`);

Create frames for an animation of the problem. (This may take some time)
for t from 0 to 60 by 1/3 do
g.(t*3):=plot([[0,0],[Qx,Qy],[Px,Py]]):
od:

Display the animation.
display([g.(0..180)], insequence=true);

Reset t variable
t:='t';

Find the velocity and speed of point P
vPx:=diff(Px,t);
vPy:=diff(Py,t);
speed:=sqrt(vPx^2+vPy^2);
speed:=simplify(speed);

plot(speed,t=0..50,title=`Speed of P vs time`);

dspeed:=simplify(diff(speed,t));

Fsolve finds a solution within the interval specified.
Time:=fsolve(dspeed=0,t=2..6);

Location of P at maximum speed.
locationP:=simplify(subs(t=Time,[Px,Py]));
locationQ:=simplify(subs(t=Time,[Qx,Qy]));

Convet to polar coordinates.

Fastest speed is when P is closest to the center.
locP:=polar(locationP[1]+locationP[2]*I);
locQ:=polar(locationQ[1]+locationQ[2]*I);

plot speed on z axis and position on x and y axises.
spacecurve([Px,Py,speed],t=0..120,numpoints=500,title=`Speed of P vs
position`);

Find the acceleration of point P
aPx:=diff(vPx,t);
aPy:=diff(vPy,t);
MagAcc:=sqrt(aPx^2+aPy^2);
MagAcc:=simplify(MagAcc);

plot(MagAcc,t=0..50,title=`|a| vs time`);

dMagAcc:=simplify(diff(MagAcc,t));

Fsolve finds a solution within the interval specified.
Time:=fsolve(dspeed=0,t=15..19);

Location of P at maximum acceleration.
locationP:=simplify(subs(t=Time,[Px,Py]));
locationQ:=simplify(subs(t=Time,[Qx,Qy]));

Acceleration is a maximum when P is furthest from the center.

Convet to polar coordinates.
locP:=polar(locationP[1]+locationP[2]*I);
locQ:=polar(locationQ[1]+locationQ[2]*I);

plot |a| on z axis and position on x and y axises.
spacecurve([Px,Py,MagAcc],t=0..120,numpoints=500,title=`|a| of P vs
position`);

PHYSICS RESOURCE PACKETS PROJECT            File: TENNIS.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A tennis player during her service throws a 59-g tennis ball straight up and at
the peak of its trajectory hits it with her racket which exerts a horizontal
force given by F = at + bt^2, where a= 1200 N/ms, and b=400 N-ms^(-2).  The ball
separates from the racket after t=3.0 ms, where t is measured from the time the
racket first contacts the ball.
a) Sketch the impulse as a function of the time using a spreadsheet and
determine i) the impulse and ii) the average impulsive force
b) Using maple or mathematica or other means, determine the impulse and the
average impulsive force.

Solution:
restart;
Enter constants
m:=59/1000; a:=1200; b:=400;

Equation for force on tennis ball.
F:=a*t-b*t^2;

plot the force on the tennis ball
plot(F,t=0..3,title=`Force vs time`);

Find the impulse on the ball.
impulse:=Int(F,t=0..3);

Evaluating the intergal can be done in two ways
Impulse:=value(impulse);
Impulse:=int(F,t=0..3);

Find the average impulse.
AvgImpulse:=Impulse/(3-0);

PHYSICS RESOURCE PACKETS PROJECT            File: WALLFALL.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A stud wall has been built by construction workers and carelessly left standing
upright unsupported. It has a mass of 11.5 kg and a height of 2.66 m. The wall
CM is 1.33 m from either end, and its moment of inertia with respect to an axis
running through the side on the ground is 36.7 kg-m^2.
The wall is accidentally brushed by a worker as he goes by, giving it an initial
a) Determine the angular velocity when the wall hits the ground.
b) Find the time for the wall to hit the ground, using a numerical integration
technique.

Solution:
restart;
Part A

Conservation of energy equation.
f:=m*g*L/2+Ip*omega[0]^2/2=Ip*omega^2/2+m*g*L*cos(theta)/2;

Define constants
m:=115/10; L:=266/100;Ip:=367/10;g:=98/10;omega[0]:=5/100;
f2:=f;

solve f2 for omega.
sol:=solve(f2,omega);

Pick positive result.
f3:=omega=sol[1];

Wall hits ground when theta=Pi/2.  Substitute theta=Pi/2 into f3.
omega_at_ground:=evalf(subs(theta=Pi/2,f3));

Part B

convert to differential form
f4:=subs(omega=d*theta/dt,f3);

solve for dt
f5:=dt=solve(f4,dt);

Remove the differentials so Maple can integrate the equation.
f6:=lhs(f5/dt)=rhs(f5/(d*theta));

Integrate both sides of the equation.  Right side of equation is integrated
numerically.
ans:=int(lhs(f6),t1=0..t)=evalf(int(rhs(f6),theta=0..Pi/2));

PHYSICS RESOURCE PACKETS PROJECT            File: YANK.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

A mass M hangs from a fixed support by a thin thread.  An identical thread is
connected to the bottom of mass M.  By grasping the bottom thread, a person can
would be done. [1-2 min TPS; yank sharply on lower thread to break it first.]

We will model the 'inertial' effect of the mass M when we pull suddenly on the
constant k, the same unstretched length Lo, and the same breaking strength B.
The threads (springs) will break when the force exceeds B and the extension of
the spring exceeds D= B/k.  D is then the maximum extension of the spring before
the string breaks.

Before we pull on the lower spring, (t<=0) mass M is at equilibrium.  The
upper spring is stretched by an amount Mg/k, and exerts an upward force of Mg,
opposite to the mass's weight. All distances are measured positive downward from
the fixed support.  Pulling on the bottom thread means that the free end  of the
lower spring travels at a constant velocity V.

h(t) is the upper spring length, from the support to mass M. At t=0,
h(0) =Lo + Mg/k.  x(t) is the distance downward from h(0) to mass M.  At t=0,
x(0)=0.
After t=0,  h(t) = Lo + Mg/k + x(t).

s(t) is the lower spring length from mass M to the free end. At t=0, s(0) = Lo.
After t=0,  s(t) = Lo+Vt -x(t).

The net force on the mass M is then Fnet = +Mg -k( h(t) - Lo ) + k( s(t) - Lo ).

a) Collect the information above, and write Fnet in terms of x(t), M, g, k, V
and t.
b) If the yank on the lower thread is violent enough, the mass will hardly move
before the lower thread snaps. This means that x(t) is practically zero, so
the force in the lower spring is only due to the rapid stretching from Vt.
Write down the net force on the mass in terms of t, assuming that x(t) is
nearly zero.
c) Solve Fnet = Ma  for the displacement x(t), assuming that in Fnet  x can be
set to zero. (The second time derivative of  x(t) won't be close to zero,
even though x is nearly zero.). [Use dsolve, or integrate by hand].

Solution:
restart;
Part A

Enter in eqautions for h(t), and s(t).
ht:=L[0]+m*g/k+xt;
st:=L[0]+V*t-xt;

Equation for Fnet.
Fnet:=m*g-k*(ht-L[0])+k*(st-L[0]);

Part B

Equation for Fnet assuming that x(t) is practically 0.
Fnet2:=subs(xt=0,Fnet);

Part C

Put equation in differential form.
FnetDE:=diff(x(t),t\$2)=Fnet2;

Solve for x(t) using initial conditions.
sol:=dsolve({FnetDE,x(0)=0,D(x)(0)=0},x(t));

Set xt equal to the "right hand side" of sol.
xt:=rhs(sol);

Problem can also be solved using integration.

Integrate once to find the velocity of the mass.  [note: maple does not add the
constant of integration, but it is not needed in this case.]
vt:=int(Fnet2,t);

Integrate a second time to find the position of the mass.  This is the same
result obtained from the differential equation.
xt:=int(vt,t);

Part D

Enter in values given in problem.
g:=98/10; m:=115/100; k:=100; L[0]:=18/100; V:=14; Fbreak:=12;

Force on lower spring.
Flower:=k*(st-L[0]);

Find the time for the lower spring to break by solving for t when Flower=Fbreak.
sol2:=fsolve(Flower=Fbreak,t);

Pick smallest posivite solution for the time that the spring will break.
TimeBreak:=sol2[2];

Part E

Equation for the position of the mass with constants.
xt;

Find the position of the mass when t=TimeBreak.

[Note that the mass has moved vary little, so our assumption is correct.]
xtBreak:=subs(t=TimeBreak,xt);

Find the force on the upper spring by multiplying k by the distance the spring
is stretched from its free length.
Fupper:=k*(xtBreak+m*g/k);

Answer can also be found by substituting in t=TimeBreak into the equation below.
Fupper:=k*(ht-L[0]);

Force on the upper string when the lower string breaks.
subs(t=TimeBreak,Fupper);

The lower string breaks first because the force on the upper string is less than
12N.

PHYSICS RESOURCE PACKETS PROJECT            File: YANKADV.MWS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

In the "yank.ms" problem, don't assume that x(t) is nearly zero in the
force equation.  Solve the full Fnet = Ma equation for x(t) using dsolve.
a) Obtain the solution x(t) in terms of two arbitrary constants.  Then use the
condition that x(0) = 0 and dx/dt(0) = 0 to evaluate the two arbitrary
constants.
b) Use the dsolve command again, but this time include the boundary conditions
right in the command, so that the solution automatically fits the two
boundary conditions.
c) Now let M = 1.15 kg,  k = 100 N/m, Lo = 0.18 m and V = 14.0 m/s. Assume
either thread will break when the force is 12.00 N.  [The upper thread is
quite close to breaking before the lower thread is pulled at all.] Calculate
the time for the lower thread to break.
d) Determine the force on the upper thread when the lower one breaks.  Which
Going Further: Solve for the velocity  V which breaks both strings at once.
Velocities faster than this should break the lower spring first, and slower
velocities will break the upper one first. [Answer: both break @ 2.13 m/s]

Solution:
restart;
with(plots):

Tell maple that k and t are positive.
assume(t>0);
assume(k>0);

Part A

Enter in eqautions for h(t), and s(t).
ht:=L[0]+m*g/k+x(t);
st:=L[0]+V*t-x(t);

Equation for Fnet.
Fnet:=m*g-k*(ht-L[0])+k*(st-L[0]);

Put equation in differential form.
FnetDE:=diff(x(t),t\$2)=Fnet;

Solve for x(t) using the dsolve command.
sol:=dsolve(FnetDE,x(t));

Set xt equal to the "right hand side" of sol.
xt:=rhs(sol);

v(t) is the derivative of x(t)
vt:=diff(xt,t);

Write two independent equations that satisfy the initial conditions of the
problem.

x(0)=0
eq1:=simplify(subs(t=0,xt)=0);

v(0)=0
eq2:=simplify(subs(t=0,vt)=0);

Solve the two equation for the two unknowns: _C1 and _C2
constants:=solve({eq1,eq2},{_C1,_C2});

Assign values to _C1 and _C2.
assign(constants);

Display the current contents of xt.
xt;

Part B

Solve for x(t) using initial conditions.
sol:=dsolve({FnetDE,x(0)=0,D(x)(0)=0},x(t));

Solution is the same as above.
Part C

Force on lower spring.
Flower:=k*(st-L[0]);

Substitute in the solution for x(t) into the above equation
Flower1:=subs(x(t)=xt,Flower);

Substitute in the values given in part C
Flowersub:=subs({k=100,V=14},Flower1);

Find the time for the lower spring to break by solving for t when
Flowersub=12 N.
TimeBreak:=fsolve(Flowersub=12,t);

Part D

Equation for the position of the mass with constants.
xtsub:=subs({V=14,k=100},xt);

Force on the upper spring.
Fupper:=k*(ht-L[0]);

Substitute in know values into the above equation.
Fuppersub:=subs({k=100,m=115/100,g=98/10,x(t)=xtsub},Fupper);

Force on the upper string when the lower string breaks.
evalf(subs(t=TimeBreak,Fuppersub));

The lower thread breaks first because the force on the upper thread is less
than 12N.

Going Further:

Find the position of the mass in terms of V.
xtsub2:=subs(k=100,xt);

The force on the upper spring in term of V is:
Fuppersub2:=subs({k=100,m=115/100,g=98/10,x(t)=xtsub2},Fupper);

Force on the lower spring in terms of V:
Flowersub2:=subs({k=100,m=115/100,g=98/10,x(t)=xtsub2},Flower);

Draw an animation of the force on the upper and lower springs as V increases.
The horizontal line on the graph is the 12N breaking point of the spring.
animate({Fuppersub2,Flowersub2,12},t=0..0.2,V=0..4,color=red,title=`Force
on the upper and lower springs vs time as speed increases`);

We want both springs to break under a 12N force at the same time.
eq1:=Fuppersub2=12;
eq2:=Flowersub2=12;

Solve the above system of equations for V and t.
sol:=fsolve({eq1,eq2},{V,t},{V=0..14,t=-1..1});

Both springs will break at the same time when V=2.13185 m/s.

Substitute in the value obtained for V into the spring force equaitons.
sub2:=subs(V=2.131848184,{Fuppersub2,Flowersub2,12});

plot(sub2,t=0..0.2,title=`Graph of forces on springs when V=2.13185 m/s`);