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

MAPLE text, input lines and file downloads (plus 1 Mma download)

bounce4.mws

A large ball and small ball fall together so that after the large ball hits the floor, it impacts the small ball and sends the small ball high into the air. This can be demonstrated on a ringstand with masses and springs. The large and small masses are separated by a spring, and there is a spring below the large mass.

This worksheet models that situation. It has a spring between the large and small mass, and a spring below the large mass, between it and the floor. It illustrates use of Heaviside functions to model contact between bodies with a spring. Motions of the bodies are found as a function of time.

The idea is to determine the separation between large and small mass for the largest rebound height of the small mass.

y1 and y2 released from rest, y1 above y2.

y1 > y2 y1 is above y2.

k=2000 N/m, m1 = 1, m2 = 2

restart;with(plots):

Spring compression. Force comes on when compression is positive

c2:=-y2(t)+L;
c1:=L+y2(t)-y1(t);
eq1:=+k*c1*Heaviside(c1)-m1*98/10 = m1*diff(y1(t),t,t);
values:={m1=1/5,m2=2,k=2000,L=1/2};
eq1s:=subs(values,eq1);
eq2:=+k*c2*Heaviside(c2)-k*c1*Heaviside(c1) -m2*98/10 = m2*diff(y2(t),t,t);
eq2s:=subs(values,eq2); > eqs:={eq1s,eq2s};
s:=dsolve({eq1s,eq2s,y1(0)=5,D(y1)(0)=0,y2(0)=37/10,D(y2)(0)=0},{y1(t),y2(t)},numeric); h:=odeplot(s,[t,y1(t)],0..7/2,color=red):
j:=odeplot(s,[t,y2(t)],0..7/2,color=blue):
display({j,h},title=`Bounce of y1(red), y2(blue) `);

Download bounce4.ms [Use 'File Save']
Download bounce4.mws [Use 'File Save']

conserve.ms

Software: Maple V Release 3.0 Authors: Perry Peters and Greg Williby
Concepts Used: Spherical coordinates

Here are the components of a supposedly conservative force

a) Convert Fx, Fy, and Fz into spherical polar coordinates using
x = r sin theta cos phi, y = r sin theta sin phi, and z = r cos theta.

b) For integrating in a circle at constant z, both r and theta must be constant.
Show that an element of length dl along a path at constant r and theta will be

This must be done by means of a sketch, showing dl, and the angles theta and phi . c) Integrate the dot product of the force F with the length element dl around a path at constant z, holding theta and r constant, showing that the value of the integral will be zero for all values of r and theta. [The integral around any closed path should be zero for a conservative force.]

Enter in the force vector.

Show that the potential energy of F is:

Perform partial derivative tests on U to see if it is the potential energy of F

U is thus the potential energy of F.
Convert the vector into polar form.

Direction of travel.

Take the dot product of the force and the direction before integrating.

Integrate around the loop.

This shows that the work done on the particle around the loop for any constant value of phi and rho is 0.

Download conserve.ms [Use 'File Save']
Download conserve.mws [Use 'File Save']
Download conserve.ma [Use 'File Save']

contact.ms

A program to illustrate use of Heaviside functions to model contact between bodies with a spring.

A moving mass m1 has a spring of unstretched length x0 attached to it. m1 has an initial velocity in the +x direction and travels until it encounters a stationary mass m2. The spring compresses during the collision until the masses separate, and we plot the motion of both masses as a function of time.

Thanks to Joaquin Casellas for repairing a flaw in the problem! (The original problem has the bodies passing through one another),

With this as a starter problem, students could be asked to work out the problem of two masses falling vertically, one above the other. One spring is attached to each mass, and one mass is directly above the other. The masses are released from rest at the same time, with some separation between them. The lower mass's spring strikes the ground and it rebounds, then encounters the spring of the upper mass. After the collisions the upper mass rebounds to some maximum height. The idea is to find the maximum rebound height of the upper mass, given masses and spring constants. The adjustable parameter is the distance between the masses on release from rest.

a) Given masses, spring constants, and unstretched spring lengths, and the initial height of the lower mass, reason out how the collision should occur so that the upper mass rebounds to a maximum height.

b) After the qualitative considerations are worked out, simulate the drop, and the motion of the two masses.

Contact between a moving mass and a stationary mass via a spring.

Download contact.ms [Use 'File Save']
Download contact.mws [Use 'File Save']

couposc.mws

coupled oscillators 7/8/96   This 'c' version will do two coupled oscillators between walls

It also teases out the normal modes (x2/x1 at each mode freq)

restart: with(plots):
L:= m1* v1^2/2 +m2*v2^2/2- k1* x1^2 / 2 -kc*(x2-x1)^2/2-k1*x2^2/2 ;
d1:=diff(L,v1 ); d2:=diff(L,v2);
p1:=subs(v1=v1(t),d1); p2:=subs(v2=v2(t),d2);

First do the derivatives, then tell it you have a function of t.

d1:=diff(p1,t)-diff(L,x1)=0;

Same thing for derivative with respect to x1. First the derivative wrt to x1, then tell it that x1 is a function of t.

d2:=diff(p2,t)-diff(L,x2)=0;
f1:=subs(x1=x1(t),x2=x2(t),v1(t)=diff(x1(t),t),d1);
f2:=subs(x1=x1(t),x2=x2(t),v2(t)=diff(x2(t),t),d2);
h1:=subs(v1(t)=diff(x1(t),t),f1);
h2:=subs(v1(t)=diff(x1(t),t),f2);
z:=lhs(h1);

Solve for the simple harmonic oscillator equation.

sol:=dsolve({h1,x1(0)=0,D(x1)(0)=1,h2,x2(0)=0,D(x2)(0)=0},{x1(t),x2(t)});
sb:={m1=1,m2=1,k1=1,kc=1/10};
sol1:=subs(sb,sol);
q1:=simplify(evalf(rhs(sol1[1])));
q2:=simplify(evalf(rhs(sol1[2])));
plot({q1,q2},t=0..50,title=`coupled oscillators`);
hx:=subs(diff(x1(t),t,t)=-omega^2*x1(t),h1);
x21:=solve(hx,x2(t));
hy:=subs(diff(x2(t),t,t)=-omega^2*x2(t),h2);
x22:=solve(hy,x2(t));
eqn2:= r*x1(t) = x22;
ratio1:=solve(eqn2,r);

This puts us in position to substitute frequencies and find the mode configurations!

solfreq:=solve(x22=x21,omega);
frq1:=solfreq[1]; > frq2:=solfreq[3];
omega1:=evalf(subs(sb,frq1));
omega2:=evalf(subs(sb,frq2));
mode2x2_over_x1:=subs(sb,omega=omega2,ratio1);

In mode 2, the masses move together, and the center spring is not compressed or expanded. The amplitude ratio is 1:1

mode1x2_over_x1:=subs(sb,omega=omega1,ratio1);

In mode 1, the masses move oppositely, compressing the spring, and operating at a higher frequency. The amplitude ratio is -1:1.

The 'transfer time' from one oscillator to the other should be argued from 'beats', where the plot looks just the same. Want 1/4 of the 'period' of the difference frequency, (f1-f2)/2

Marion problem 12-5, different masses, same k

solm:=solve((-m*x+2*k)*(-M*x+2*k) = k^2,x);
A_B_ratio:=(-m*x+2*k)/k;
rat1:=subs(x=solm[1],A_B_ratio);
rs:=simplify(rat1);
subs(M=10*m,rs);

Download cuploscc.ms [Use 'File Save']
Download cuploscc.mws [Use 'File Save']

drivnosc.ms

Plot the phase and amplitude response of a driven, damped oscillator as a function of frequency.
It is interesting that very heavy damping (b=10, k=10, m=1) crushes the response away from zero frequency.

Phase of the response as a function of frequency. The conventional definition is the phase of the driver minus the phase of the driven system. The relative amplitude response will be taken as exp(-I alph).

The driver is in phase with the driven system at very low frequencies, and pulls ahead until it is Pi/2 ahead at resonance, and goes on to be Pi ahead at very high frequencies.

Amplitude response as a function of frequency

Download drivnosc.ms [Use 'File Save']
Download drivnosc.mws [Use 'File Save']

eulerAn2.mws           6/23/97

Show the effect of a rotation matrix on a vector as a function of time. In particular, implement the euler angle rotations and animate them.

Viewing is done along a polar (Pol) and azimuthal (Az) angle.

Initial angles theta, phi, and psi are specified, and rotation rates for all 3.

The calculation starts from body x3,y3,z3 components, and works backward to the space xyz coordinates. When x,y,z are known for each body axis, a projection is done into the viewing direction, and an animation is made.

Forward rotation thru angle a gives xnew = xold cos(a) + yold sin(a), etc.

restart;with(plots):
with(linalg):

m1 is forward rotation thru x. new coords = m1 * old coords ditto for m2 and m3

m1:=matrix([[1,0,0],[0,cos(x),sin(x)],[0,-sin(x),cos(x)]]);
m2:=matrix([[cos(x),0,sin(x)],[0,1,0],[-sin(x),0,cos(x)]]);
m3:=matrix([[cos(x),sin(x),0],[-sin(x),cos(x),0],[0,0,1]]);

We start with body'x axes in x3,y3,z3 system as (1,0,0) for x axis, etc. To get to components in original space xyz system we must do a reverse rotation thru -psi, then reverse thru -theta, then reverse thru -phi. This will lead us to the xyz components of the body's axes after the body has been rotated in space

r1:=subs(x=-phi,evalm(m3));
r2:=subs(x=-theta,evalm(m1));
r3:=subs(x=-psi,evalm(m3));
r123:=evalm(r1 &* r2 &* r3);
v1:=matrix([[f1],[f2],[f3]]);
v:=evalm(r123 &* v1);
p0:=[0,0]; > pL:=[-a*sin(Az)+b*cos(Az),-cos(Pol)*(a*cos(Az)+b*sin(Az))+c*sin(Pol)];

Viewing angles are set up here, polar and azimuthal

viewsubs:={Pol=2*Pi/10,Az=1*Pi/4};

Now set up the view and plots of the x,y,z axes

ax:=subs(a=1,b=0,c=0,viewsubs,pL):
ay:=subs(a=0,b=1,c=0,viewsubs,pL):
az:=subs(a=0,b=0,c=1,viewsubs,pL):
zax:=plot([p0,az],color=black,axes=NONE,scaling=constrained):
xax:=plot([p0,ax],color=black,axes=NONE,scaling=constrained): yax:=plot([p0,ay],color=black,axes=NONE,scaling=constrained):

Now set up the view of lines which have been rotated

p1:=subs(a=v[1,1],b=v[2,1],c=v[3,1],pL):
p:=subs(viewsubs,p1): > px:=subs(f1=1,f2=0,f3=0,p):
py:=subs(f1=0,f2=1,f3=0,p): > pz:=subs(f1=0,f2=0,f3=1,p):
Nframes:=35; > Psidot:=2*Pi/18; Thetadot:=0;Phidot:=2*Pi/36;
Theta0:=Pi/12;Phi0:=0;Psi0:=0;

for i from 0 to Nframes do
sbs:={theta=Theta0+i*Thetadot,phi=Phi0+i*Phidot,psi=Psi0+i*Psidot}: gx.(i):=plot(subs(sbs,[p0,px]),color=blue,scaling=constrained,thickness=3,axes=NONE):
gy.(i):=plot(subs(sbs,[p0,py]),color=red,scaling=constrained,thickness=3,axes=NONE):
gz.(i):=plot(subs(sbs,[p0,pz]),thickness=3,color=black,scaling=constrained,axes=NONE):
m.(i):=display({gx.(i),gy.(i),gz.(i),xax,yax,zax}): od:
display([m.(0..Nframes)],insequence=true);

Download eulerAn2.mws [Use 'File Save']

fouclt.mws

Foucault pendulum simulation - nominal values Consult Marion/Thornton circa p. 400 for treatment

restart;with(plots):
Rotation rate is set at 1/10 to make things visible on a short plot

omega:=1/10;
eq1:=diff(x(t),t,t)+98/20*x(t) = omega*diff(y(t),t);
eq2:=diff(y(t),t,t)+(98/20)*y(t) = -omega*diff(x(t),t); s:=dsolve({eq1,eq2,x(0)=1/10,D(x)(0)=0,y(0)=0,D(y)(0)=0},{x(t),y(t)},numeric);
xp:= z -> rhs(s(z)[2]);
yp:= z -> rhs(s(z)[4]);
plot([xp,yp,0..10],title=`x-y plot of motion`);

Download fouclt.mws [Use 'File Save']

The following file is a template for students, needing them to fill in equations to get the plot
Download foucmt.mws [Use 'File Save']

proj.ms

Projectile motion in Lagrangian formalism

.
An introductory experience with lagrangian mechanics in a familiar setting.


Valuable as a template for making Maple do what you want it to.

Download proj.ms [Use 'File Save']
Download proj.mws [Use 'File Save']

rotspr.ms

Two equal masses are connected by a spring, and are rotating about their CM.

rotating mass-spring system in Lagrangian formalism, patterned after proj.ms Maple seems to need a lot of shifting around in the variables.

Generalized momenta w/ symbols (pj = dLdqjdot)

The square brackets are important in the following dsolve command to keep the solutions in correct order

Download rotspr.ms [Use 'File Save']
Download rotspr.mws [Use 'File Save']

table.ms

A mass m on a frictionless tabletop is connected by a string of length H through a hole in the table to a mass M hanging below the table.

Done via lagrange's equations

Students might be asked to

a) Identify the two constants of the motion.
b) Given m and M and initial r and phi dot values, predict the maximum r analytically.
c) Check the answer with the maple solution.
d) Determine the maximum tension in the cord from the conditions of the problem.

The square brackets are important in the following dsolve command to keep the solutions in correct order

Download table.ms [Use 'File Save']
Download table.mws [Use 'File Save']

main line problems