Application Center - Maplesoft

# Motion of a Heavy Symmetric Top

You can switch back to the summary page by clicking here.

heavytop.mws

Motion of a heavy symmetric top

Frank Wang and J. F. Ogilvie

The Cooper Union for the Advancement of Science and Art

Department of Physics, Cooper Square, New York, NY 10003

wang3@cooper.edu

Centre for Experimental and Constructive Mathematics, Department of Mathematics, Simon Fraser University,
8888 University Drive, Burnaby, British Columbia V5A 1S6 Canada

Calculus of variations involves finding a derivative of a function with respect to a function, which Maple does not directly support. With simple substitution, we can, however, solve easily a problem of this type without invoking an external library. We use a problem in classical mechanics to illustrate a procedure to solve the Euler-Lagrange equation, which is the most important application of calculus of variations.

A symmetric top with one point fixed is under the influence of gravity. Let the symmetry axis be taken as z fixed in the body; the moment of inertia is I 3 along this axis. By symmetry I 1 = I 2. The configuration of the top is completely specified by three Euler's angles: gives the inclination of axis z from the vertical, measures the azimuth of the top about the vertical, and is the angle of rotation of the top about its own z axis. Please refer to Goldstein's Classical Mechanics (1980), p. 213, for definition and detailed derivation. The lagrangian, which is kinetic energy minus potential energy, can be expressed in terms of those eulerian angles.

We find the equations of motion by employing the Euler-Lagrange equation:

- = 0

in which is a coordinate and denotes its derivative with respect to time: = ; conventional notation for this quantity as with a dot above this letter is impracticable in Maple text.

> restart:

> T := I1/2*(diff(theta(t),t)^2 + diff(phi(t),t)^2*sin(theta(t))^2)
+ I3/2*(diff(psi(t),t) + diff(phi(t),t)*cos(theta(t)))^2;

> V := m*g*l*cos(theta(t));

The lagrangian is kinetic energy minus potential energy.

> L := T - V;

Six variables in this lagrangian -- , , and their derivatives -- are functions of time that can not be differentiated directly in Maple ; we therefore substitute them with symbols var1 , var2 etc.

> L1 := subs({phi(t)=var1, theta(t)=var2, psi(t)=var3, diff(phi(t), t)=var4,
diff(theta(t), t)=var5, diff(psi(t), t)=var6}, L);

This lagrangian contains no explicit dependence on and ; the corresponding generalized momenta are constant in time.

For coordinate , the momentum is

> Eq11 := diff(L1, var6);

We can substitute var1 , var2 ... that we defined earlier back to their original assignments.

> Eq13 := subs({var1=phi(t), var2=theta(t), var3=psi(t), var4=diff(phi(t), t),
var5=diff(theta(t), t), var6=diff(psi(t), t)}, Eq11);

We define this constant momentum as Mz .

> Eq14 := Eq13 = Mz;

We proceed similarly for coordinate ; we define this constant momentum to be Mz1 .

> Eq21 := diff(L1, var4);

> Eq23 := subs({var1=phi(t), var2=theta(t), var3=psi(t), var4=diff(phi(t), t),
var5=diff(theta(t), t), var6=diff(psi(t), t)}, Eq21);

> Eq24 := Eq23 = Mz1;

Applying the Euler-Lagrange equation to coordinate r , we obtain

> Eq31 := diff(L1, var5);

> Eq32 := diff(L1, var2);

> Eq33 := subs({var1=phi(t), var2=theta(t), var3=psi(t), var4=diff(phi(t), t),
var5=diff(theta(t), t), var6=diff(psi(t), t)}, Eq31);

> Eq34 := subs({var1=phi(t), var2=theta(t), var3=psi(t), var4=diff(phi(t), t),
var5=diff(theta(t), t), var6=diff(psi(t), t)}, Eq32);

> Eq35 := diff(Eq33, t);

> Eq36 := Eq35 - Eq34 = 0;

We list three equations of motion.

> Eq14;

> Eq24;

> Eq36;

To these equations there is no analytic solution, but when we supply initial conditions we can solve these equations numerically. It is convenient to express these equations in terms of constant Mz and Mz1 .

From equations Eq14 and Eq24 we eliminate to find .

> Eq51 := isolate(Eq14, diff(psi(t), t));

> Eq52 := isolate(Eq24, diff(psi(t), t));

> Eq53 := rhs(Eq51) = rhs(Eq52);

> Eq54 := isolate(Eq53, diff(phi(t), t));

Substituting back into either Eq14 or Eq24 , we find .

> Eq55 := subs(Eq54, Eq51);

Finally substituting and into Eq36 , we obtain an equation containing a single variable .

> Eq61 := subs({Eq54, Eq55}, Eq36);

We list three equations in terms of .

> Eq54;

> Eq55;

> Eq61;

For a numerical worked example we assume a top constructed from a heavy circular disk of mass = 0.1 kg and radius = 0.02 m mounted at the center of a thin rod of length 0.02 m. Suppose that we start the top spinning at 35 rev s , and place it at angle rad = 60 deg from the vertical. At we assume = 0 and = -5 rad s .

> g := 9.8;

> m := 0.1;

> l := 0.02;

> R := 0.02;

The moment of inertia of a disk is and .

> I1 := m*l^2;

> I3 := 1/2*m*R^2;

is given, along with other inital conditions.

> Eq71 := diff(psi(t), t) = evalf(35*2*Pi);

> Eq72 := phi(0) = 0;

> Eq73 := D(phi)(0) = -5; #try different values

> Eq74 := theta(0) = evalf(Pi/3);

> Eq75 := D(theta)(0) = 0;

> Mz := subs({Eq71, diff(phi(t), t) = rhs(Eq72)}, lhs(Eq14));

> Mz1 := evalf(subs({Eq71, diff(phi(t), t)=rhs(Eq73), theta(t)=rhs(Eq74)}, lhs(Eq24)));

Here are these equations again with supplied initial conditions.

> Eq54;

> Eq61;

We solve Eq54 and Eq61 numerically, with initial conditions Eq72 , Eq74 and Eq75 .

> Eq91 := dsolve({Eq54, Eq61, Eq72, Eq74, Eq75}, {theta(t), phi(t)},
numeric, output=listprocedure);

We plot the shapes for the locus of the figure axis on a unit sphere.

> with(plots):

Warning, the name changecoords has been redefined

> p1 := spacecurve([l, rhs(Eq91(t)[2]), rhs(Eq91(t)[3])], t=0..0.5,
coords=spherical, numpoints=100, thickness=2, colour=black):

> p2 := sphereplot(l, theta=0..2*Pi, phi=0..Pi, style=hidden):

> display([p1,p2], scaling=constrained);

For another condition we take at .

> Eq73 := D(phi)(0) = 0;

> Mz;

> Mz1 := evalf(subs({Eq71,diff(phi(t), t) = rhs(Eq73), theta(t) = rhs(Eq74)}, lhs(Eq24)));

> Eq54:

> Eq61:

> Eq91 := dsolve({Eq54,Eq61,Eq72,Eq74,Eq75}, {theta(t),phi(t)}, numeric,
output=listprocedure):

> p3 := spacecurve([l, rhs(Eq91(t)[2]), rhs(Eq91(t)[3])], t=0..0.5,
coords=spherical, numpoints=100, thickness=2, colour=black):

> display([p3, p2], scaling=constrained);

For a further initial condition, we take at .

> Eq73 := D(phi)(0) = 2;

> Mz;

> Mz1 := evalf(subs({Eq71, diff(phi(t), t) = rhs(Eq73), theta(t) = rhs(Eq74)}, lhs(Eq24)));

> Eq54:

> Eq61:

> Eq91 := dsolve({Eq54, Eq61, Eq72, Eq74, Eq75}, {theta(t), phi(t)},
numeric, output=listprocedure):

> p4 := spacecurve([l, rhs(Eq91(t)[2]), rhs(Eq91(t)[3])], t=0..0.5,
coords=spherical, numpoints=100, thickness=2, colour=black):

> display([p4, p2], scaling=constrained);

>