Application Center - Maplesoft

App Preview:

Motion of a Heavy Symmetric Top

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

Learn about Maple
Download Application


 

heavytop.mws

Motion of a heavy symmetric top

Frank Wang ``^`1` and J. F. Ogilvie ``^`2`

``^`1` The Cooper Union for the Advancement of Science and Art

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

wang3@cooper.edu

``^`2` 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: theta gives the inclination of axis z from the vertical, phi measures the azimuth of the top about the vertical, and psi 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:

d/dt diff(L,`q'`[i]) - diff(L,q[i]) = 0

in which q is a coordinate and q*`'` denotes its derivative with respect to time: q*`'` = dq/dt ; conventional notation for this quantity as q 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;

T := 1/2*I1*(diff(theta(t),t)^2+diff(phi(t),t)^2*si...

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

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

The lagrangian is kinetic energy minus potential energy.

> L := T - V;

L := 1/2*I1*(diff(theta(t),t)^2+diff(phi(t),t)^2*si...

Six variables in this lagrangian -- theta , phi , psi 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);

L1 := 1/2*I1*(var5^2+var4^2*sin(var2)^2)+1/2*I3*(va...

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

For coordinate psi , the momentum is

> Eq11 := diff(L1, var6);

Eq11 := I3*(var6+var4*cos(var2))

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);

Eq13 := I3*(diff(psi(t),t)+diff(phi(t),t)*cos(theta...

We define this constant momentum as Mz .

> Eq14 := Eq13 = Mz;

Eq14 := I3*(diff(psi(t),t)+diff(phi(t),t)*cos(theta...

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

> Eq21 := diff(L1, var4);

Eq21 := I1*var4*sin(var2)^2+I3*(var6+var4*cos(var2)...

> 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);

Eq23 := I1*diff(phi(t),t)*sin(theta(t))^2+I3*(diff(...

> Eq24 := Eq23 = Mz1;

Eq24 := I1*diff(phi(t),t)*sin(theta(t))^2+I3*(diff(...

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

> Eq31 := diff(L1, var5);

Eq31 := I1*var5

> Eq32 := diff(L1, var2);

Eq32 := I1*var4^2*sin(var2)*cos(var2)-I3*(var6+var4...

> 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);

Eq33 := I1*diff(theta(t),t)

> 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);

Eq34 := I1*diff(phi(t),t)^2*sin(theta(t))*cos(theta...

> Eq35 := diff(Eq33, t);

Eq35 := I1*diff(theta(t),`$`(t,2))

> Eq36 := Eq35 - Eq34 = 0;

Eq36 := I1*diff(theta(t),`$`(t,2))-I1*diff(phi(t),t...

We list three equations of motion.

> Eq14;

I3*(diff(psi(t),t)+diff(phi(t),t)*cos(theta(t))) = ...

> Eq24;

I1*diff(phi(t),t)*sin(theta(t))^2+I3*(diff(psi(t),t...

> Eq36;

I1*diff(theta(t),`$`(t,2))-I1*diff(phi(t),t)^2*sin(...

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 psi to find phi .

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

Eq51 := diff(psi(t),t) = Mz/I3-diff(phi(t),t)*cos(t...

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

Eq52 := diff(psi(t),t) = (Mz1-I1*diff(phi(t),t)*sin...

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

Eq53 := Mz/I3-diff(phi(t),t)*cos(theta(t)) = (Mz1-I...

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

Eq54 := diff(phi(t),t) = -(Mz*cos(theta(t))-Mz1)/I1...

Substituting phi back into either Eq14 or Eq24 , we find psi .

> Eq55 := subs(Eq54, Eq51);

Eq55 := diff(psi(t),t) = Mz/I3+(Mz*cos(theta(t))-Mz...

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

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

Eq61 := I1*diff(theta(t),`$`(t,2))-1/I1*(Mz*cos(the...

We list three equations in terms of theta(t) .

> Eq54;

diff(phi(t),t) = -(Mz*cos(theta(t))-Mz1)/I1/sin(the...

> Eq55;

diff(psi(t),t) = Mz/I3+(Mz*cos(theta(t))-Mz1)/I1/si...

> Eq61;

I1*diff(theta(t),`$`(t,2))-1/I1*(Mz*cos(theta(t))-M...

For a numerical worked example we assume a top constructed from a heavy circular disk of mass m = 0.1 kg and radius R = 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 ``^(-1) , and place it at angle theta(0) = Pi/3 rad = 60 deg from the vertical. At t = 0 we assume theta*`'` = 0 and phi*`'` = -5 rad s ``^(-1) .

> g := 9.8;

g := 9.8

> m := 0.1;

m := .1

> l := 0.02;

l := .2e-1

> R := 0.02;

R := .2e-1

The moment of inertia of a disk is m*R^2/2 and I3 = m*l^2 .

> I1 := m*l^2;

I1 := .4e-4

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

I3 := .2000000000e-4

psi*`'` is given, along with other inital conditions.

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

Eq71 := diff(psi(t),t) = 219.9114858

> Eq72 := phi(0) = 0;

Eq72 := phi(0) = 0

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

Eq73 := D(phi)(0) = -5

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

Eq74 := theta(0) = 1.047197551

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

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

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

Mz := .4398229716e-2

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

Mz1 := .2024114858e-2

Here are these equations again with supplied initial conditions.

> Eq54;

diff(phi(t),t) = -25000.00000*(.4398229716e-2*cos(t...

> Eq61;

.4e-4*diff(theta(t),`$`(t,2))-25000.00000*(.4398229...
.4e-4*diff(theta(t),`$`(t,2))-25000.00000*(.4398229...

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);

Eq91 := [t = proc (t) option `Copyright (c) 1993 by...

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);

[Maple Plot]

For another condition we take phi*`'` = 0 at t = 0 .

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

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

> Mz;

.4398229716e-2

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

Mz1 := .2199114859e-2

> 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);

[Maple Plot]

For a further initial condition, we take phi*`'` = 2 at t = 0 .

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

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

> Mz;

.4398229716e-2

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

Mz1 := .2269114858e-2

> 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);

[Maple Plot]

>