Mechanics (comments to: moloney@nextwork.rose-hulman.edu)

  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!

Download airesist.mws maple file (Use 'Save File')



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. Download aveacc.mws (use 'Save File')
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. Download avevel.mws (Use 'Save File')
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); Download ballup.mws (Use 'Save File')
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; Download conshort.mws (Use 'Save File')
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. lengthR:=1/2; lengthL:=1; radvel:=2*Pi; positionM:=rhs(subs({R=lengthR,L=lengthL,omega=radvel},f5)); x and y coordinates of wheel. Rx:=lengthR*cos(radvel*t); Ry:=lengthR*sin(radvel*t); Rx2:=lengthR*cos(radvel*t+2*Pi/3); Ry2:=lengthR*sin(radvel*t+2*Pi/3); Rx3:=lengthR*cos(radvel*t+4*Pi/3); Ry3:=lengthR*sin(radvel*t+4*Pi/3); 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); Download crank.mws (Use'Save File')
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); Download hammer.mws (Use 'Save File')
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. Start with the y-components. 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); Download jim&phil.mws (Use 'Save File')
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. Download launch.mws (Use 'Save File')
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. motion30:=combine(subs(theta=30*Pi/180,xt),radical); 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`); Download massinc.mws (Use 'Save File')
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`); Download massinfr.mws (Use 'Save File')
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. Download mspring.mws (Use 'Save File')
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)`); Download rcspring.mws (Use 'Save File')
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; Download rodgrav.mws (Use 'Save File')
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. Download skboard.mws (Use 'Save File')
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): readlib(polar): 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): readlib(polar): 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`); Download teacup.mws (Use 'Save File')
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); Download tennis.mws (Use 'Save File')
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 angular velocity of 0.05 radians/sec. 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)); Download wallfall.mws (Use 'Save File')
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 make either the top thread bread or the bottom thread break. Discuss how this 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 lower thread by treating the threads as springs. Each thread has the same spring 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. Download yank.mws (Use 'Save File')
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 thread broke first? 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`); Download yankadv.mws (Use 'Save File')

Back to Mechanics Resource Page