Unit 9: Calculus of Variations

Chapter 48: Variational Mechanics

Section 48.4: the spherical pendulum

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> libname:="C:/Program Files/Maple 6.01/Calcvar",libname:

> with(Calcvar):

> libname:="C:/Program Files/Maple 6.01/Partials",libname:
with(Partials):

> with(linalg):
with(plots):
with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

The Model

A pendulum not constrained to move in a vertical plane becomes a spherical pendulum. The pendulum bob is free to traverse the surface of a sphere centered at the pendulum's support. In fact, the spherical pendulum is the simplest pendulum to construct, since it requires only a piece of string, a bob, and a support. There is no need for a mechanism to constrain the oscillations to a single vertical plane.

If the spherical pendulum is started with motion in a single vertical plane, subsequent oscillations will remain in that plane. Hence, three-dimensional motion requires an initial velocity with a component not in the plane formed by the equilibrium position, the support, and the raised bob. As before, we will ignore friction.

Except for its moment of inertia, a small sphere such as a marble or ball-bearing rolling in a spherical bowl would be another example of the spherical pendulum.

Let us put our pendulum support at the origin, and describe the location of the bob in spherical coordinates. If the length of the pendulum is , then the coordinates of the bob are

where in measures rotation counterclockwise around the -axis, starting from the positive -axis, and in measures the angle from the positive -axis to the radius vector. (See Figure 22.8, and equations (22.10) in Section 22.3.)

>

The Lagrangian

The coordinates of the bob are given by

> xt := r*cos(u(t))*sin(v(t));
yt := r*sin(u(t))*sin(v(t));
zt := r*cos(v(t));

>

The radius vector from the origin (where the support is located) to the bob is

> R := vector([xt,yt,zt]);

>

The velocity vector for the bob is

> `R'` := map(diff,R,t);

>

A bob of mass then has kinetic energy , or

> T := collect(simplify(m/2*innerprod(`R'`,`R'`)),[m,r]);

>

Ordinary motions of the spherical pendulum find the bob below the -plane where and . Hence, for such motions, , so taking the potential energy as

> V := m*g*zt;

>

means the lowest potential is when the bob is at equilibrium, that is, when the bob is straight down. Raising the bob increases the potential energy since the potential goes to zero when the bob is lifted to the -plane. Because of the simplicity of the expression for the potential, we accept this definition, and write the Lagrangian as

> L := T-V;

>

The Equations of Motion

Using Maple's Euler_Lagrange command, we obtain

> q := Euler_Lagrange(L,t,[u(t),v(t)]);

>

We have obtained both equations of motion, namely

and

In addition, we have two first integrals. One will be the constancy of the total energy, and one is the integrated form of the first Euler-Lagrange equation. If we remove both first integrals, we are left with the two equations

> q1 := remove(has,q,[K[1],K[2]]);

>

The absence and presence of , the gravitational constant, distinguishes one equation from the other, so we single out the two equations by

> Q2 := op(select(has,q1,g)) = 0;
Q1 := op(remove(has,q1,g)) = 0;

>

To see which equation is the equation, and which is the equation, we construct the Euler-Lagrange equations from first principles using the pdiff command. For example, the term in the equation is

> pdiff(L,u);

>

because coes not explicitly contain the variable . Consequently, the equation is just , and the expanded form of this equation is

> pdiff(L,u(t)) - diff(pdiff(L,diff(u(t),t)),t) = 0;

>

Hence, the equation without is the equation, and that is the equation which has the first integral .

The equation is then

> pdiff(L,v(t)) - diff(pdiff(L,diff(v(t),t)),t) = 0;

>

and it contains the constant .

Maple's version of the first integral is

> pdiff(L,diff(u(t),t)) = c;

>

which is not the simplest form possible. Maple is programmed to prefer to the simpler , and any effort to coerce the latter form is generally wasted because, at the first opportunity, Maple will revert to the former. However, although writing

would lead to the elimination of in the equation, and hence to a quadrature in terms of elliptic functions, we take the more pragmatic approach of generating a numeric solution.

>

Example 48.7

Let us take a time-scale for which , and set both and . The Euler-Lagrange equations become

> Q3 := subs(r=1,m=1,g=1,Q1);
Q4 := subs(r=1,m=1,g=1,Q2);

>

We can obtain a numeric solution via the syntax

> QQ := dsolve({Q3,Q4, u(0)=0,v(0)=3*Pi/4, D(u)(0)=1/2,D(v)(0)=0}, {u(t),v(t)},numeric);

>

We chose initial conditions corresponding to the bob being raised through an angle of in the -plane, then being given an initial velocity in the -direction.

Maple has merely written a procedure, which, when invoked, generates the desired numeric solution. To extract this numeric solution, we write the following functions which give the values for and .

> U := s -> subs(QQ(s),u(t));
V := s -> subs(QQ(s),v(t));

>

Figure 48.11 contains graphs of in black and in cyan. It indicates that the angle is an increasing function of time, . Thus, the bob tends to rotate in a counterclockwise direction. The height of the bob is governed by the angle , and this angle remains bounded away from zero. Hence, the bob does not descend to equilibrium, but rather, oscillates between maximum and minimum heights.

> plot([U,V],0..15,color=[black,cyan],labels=[t,``], xtickmarks=7, ytickmarks=8, labelfont=[TIMES,ITALIC,12]);

>

Figure 48.12, another way of observing the motion, is a parametric plot of against . It shows that , represented by the horizontal axis, is increasing, whereas , represented by the vertical axis, undergoes a bounded oscillation with a nonzero lower bound.

> plot([U,V,0..15],color=black, labels=[u,`v `], view=[0..15,0..3], xtickmarks=7, ytickmarks=3, labelfont=[TIMES,ITALIC,12], scaling=constrained);

>

The actual motion of the bob takes place in three dimensions. Hence, Figures 48.13 and 48.14 show the motion as a space curve. In Figure 48.13, the space curve is viewed from directly above, so the curve seen is the physical trajectory's projection onto the -plane. The initial point of the trajectory is near the center of the right edge of the graph, at . The figure reveals a looping motion, with the long asis of the loops precessing counterclockwise about the -axis.

> spacecurve([cos('U'(t))*sin('V'(t)),sin('U'(t))*sin('V'(t)),cos('V'(t))],t=0..50, color=black, axes=boxed, numpoints=200, orientation=[-90,0], labels=[x,y,z],labelfont=[TIMES,ITALIC,12],tickmarks=[7,7,0]);

>

The single quote marks in the spacecurve command are required because U and V represent procedures.

Figure 48.14 shows the trajectory of the bob in space.

> spacecurve([cos('U'(t))*sin('V'(t)),sin('U'(t))*sin('V'(t)),cos('V'(t))],t=0..30, color=black, axes=boxed, numpoints=200, scaling=constrained, labels=[`x `,` y`,`z `], labelfont=[TIMES,ITALIC,12], tickmarks=[7,7,[-1,-.7]], view=[-.71..0.71,-.71..0.71,-1..-.6], orientation=[-45,45]);

>