Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Archive for March, 2010

using nm to track down the origin of an undefined symbol.

Posted by peeterjoot on March 22, 2010

The shared library for our product for historical reasons is built with the (evil) -Bsymbolic flag.

This has a number of unfortunate side effects, one of which is deferring link errors due to unresolved symbols to executable link time, instead of when we build the shared library itself.

Having screwed up in my own build over the weekend, I see a link error in my build summary this morning:

/view/peeterj_m20/vbs/engn/lib/ undefined reference to `XXYY'

I happen to know exactly what caused this (this time), but just the other day I was asked by somebody to help them figure out something similar. So the topic seems worthy of a quick blog post.

Suppose one collects all the archive names that contribute to the shared library in question (for DB2 builders one can relink the library having done an ‘export VERY_VERBOSE=true’ on the command line).

A good way to get this is to redirect that build output to a file, grab just the link line and run something like:

# cat filewithRawLinkLine.txt | tr ' ' '\n' | grep -F -e .a -e .o > archivesAndObjectNames.txt

In my case, this looked like:

# head -5 archivesAndObjectNames.txt    

This now leaves you in the position to search for the archive(s) that contributed the undefined symbol. This can be done with:

# for i in `cat archivesAndObjectNames.txt` ; do echo $i ; nm $i | grep XXYY ; done | tee nm.out | grep -B1 XXYY
 /view/peeterj_m20/vbs/engn/sqe/alibsqe.a        U XXYY

So, the archive supplying a reference to this symbol is ‘alibsqe.a’ This completes the task of doing a first localization of where the link error is coming from. One can continue to do this in a brute force way finding the object file within the archive, or solve it by source inspection once one knows where to look a bit better. In an extremely large source base, where nobody really knows any sizeable portion of it, narrowing down the problem is often an important first step.

Posted in Development environment | Tagged: , , , | 1 Comment »

Cheating header file dependencies with touch and checkin -ptime.

Posted by peeterjoot on March 19, 2010

In our build (DB2) we have some header files, that if changed will result in thousands of source files being recompiled. Sometimes you make a change to a header that you know doesn’t really require rebuilding everything. One way of cheating the build dependency system is to judiciously use the touch command to force the date backwards, as in

touch -t197301010000 /vbs/engn/include/sqleSmartArrayDefines.h
cleartool checkin -nc -ptime /vbs/engn/include/sqleSmartArrayDefines.h

An ls -l command after this will show the file time ‘1973-01-01’ (somewhat arbitrarily picked). Thirty five years or so should be early enough to not have make notice it;)

Since we use clearcase for our version control, a checkin of the file into your working branch will change the timestamp unless the -ptime option is used. Care has to be required if you are going to cheat, but I’ve seen people turning off all dependency checking in their build to try to avoid full builds after touching headers, so selecting cheating in a controlled fashion if you are going to do it is a much better option.

Another handy way to do this sort of cheating is ‘touch –reference’. This is a gnu-touch specific (ie: works on Linux) option, that allows you to specify the timestamp using a file instead of a hardcoded date, so you can match the file modification time to some previous point.


touch --reference sqle_ca_smart_array_cmd.h@@/main/temp_peeterj_db2_galileo_defect_saDeleteDebug/2 sqle_ca_smart_array_cmd.h

In this example, the checked out version (from clearcase) was touched to the predecessor version timestamp after making a change to just a comment that I didn’t want to rebuild everything for.

Posted in perl and general scripting hackery | Tagged: , , , | Leave a Comment »

stripping color control characters from lssam output.

Posted by peeterjoot on March 17, 2010

There’s probably a billion ways to do this, but here’s one that appears to work.  If you have TSA command output that has been collected by a somebody who did not use the –nocolor option.

perl -p -e 's/\x1b\[\d+m//g;' < lssamBEFORE.out

The \x1b is the ESC character.  This says to remove that ESC followed by a ‘[‘ character, then 1 or more digits and the character ‘m’, and do it for all lines.

Posted in perl and general scripting hackery | Tagged: , , , | Leave a Comment »

A matrix algebra attack on the multiple mass spherical pendulum problem. (revised)

Posted by peeterjoot on March 17, 2010

In a previous post, I’d posted the abstract and intro for a multiple mass spherical pendulum Lagrangian evaluation, along with the full pdf version of those notes.

Having noticed some typesetting bugs (braces missing after the conversion from the original Geometric Algebra version on the same topic), I’ve now posted an update. Again, it is only in pdf, since my latex to wordpress script still needs to be taught how to deal with numbering ‘subequation’s. Now also included is a summary of the equations since it was hard to see what the results and definitions were after the mess of all the algebra.

Posted in Math and Physics Learning. | Tagged: , , , , , | Leave a Comment »

vim: filtering a function body through a script

Posted by peeterjoot on March 9, 2010

While working in a coworker’s cubicle today, I casually ran the following from my vim command prompt

:,/^}/!pdcode -i 1 -s 1

I got a “whoa, what the hell was that” in response, since it was apparently not obvious. There’s really three parts to this. The first is a range expression. In vim (or vi) one can run a command on some set of lines. For example, to delete lines 4 to 8 one would do:

:4,8 d

or to delete 4 lines starting from the current position you could do (without knowing your line number)

:,+4 d

The comma by itself essentially means the current line number, so this is a starting line number of the current spot continuing for 4 additional lines down.

Now, part of the power of the line selection is the fact that one can use search expressions instead of numbers. So, I can run something like:

:,/^}/ d

and this will delete all the lines until the first line that starts with a ending curly brace. The ^ (carat) symbol is a special regular expression indicating the begining of the line, so /^}/ would not match a line that happened to have a brace somewhere else as in the following example.

if ( blah ) { }

Next ingredient is the ! (exclaimation). This says to take the selected lines and pipe them through a filter (something that takes stdin and produces modified stdout, such as grep). In my case, the script was ‘pdcode -i 1 -s 1’. This happens to be a little script I wrote, and it takes all lines like

pdTraceCodePath( 10 ) ;

in a function and renumbers them using the desired step and increment size. For reference this pdTraceCodePath is an internal db2 trace macro that we have that just adds a bit in a 64-bit word, so that later if we want a tracer of the codepath taken (all the branches) we can just dump out this single bitmask and know what the program flow was.

Posted in perl and general scripting hackery | Tagged: , | Leave a Comment »

Newton’s method for intersection of curves in a plane.

Posted by peeterjoot on March 7, 2010

[Click here for a PDF of this post with nicer formatting]Note that this PDF file is formatted in a wide-for-screen layout that is probably not good for printing.


Reading the blog post Problem solving, artificial intelligence and computational linear algebra some variations of Newton’s method for finding local minumums and maximums are given.

While I’d seen the Hessian matrix eons ago in the context of back propagation feedback methods, Newton’s method itself I remember as a first order root finding method. Here I refresh my memory what that simpler Newton’s method was about, and build on that slightly to find the form of the solution for the intersection of an arbitrarily oriented line with a curve, and finally the problem of refining an approximation for the intersection of two curves using the same technique.

Root finding as the intersection with a horizontal.

Refining an approximate horizontal intersection.

The essence of Newton’s method for finding roots is following the tangent from the point of first guess down to the line that one wants to intersect with the curve. This is illustrated in figure (\ref{fig:newtonsIntersectionHorizontal}).

Algebraically, the problem is that of finding the point x_1, which is given by the tangent

\begin{aligned}\frac{f(x_0) - b}{x_0 - x_1} = f'(x_0).\end{aligned} \hspace{\stretch{1}}(2.1)

Rearranging and solving for x_1, we have

\begin{aligned}x_1 = x_0 - \frac{f(x_0) - b}{f'(x_0)}\end{aligned} \hspace{\stretch{1}}(2.2)

If one presumes convergence, something not guarenteed, then a first guess, if good enough, will get closer and closer to the target with each iteration. If this first guess is far from the target, following the tangent line could ping pong you to some other part of the curve, and it is possible not to find the root, or to find some other one.

Intersection with a line.

Refining an approximation for the intersection with an arbitrarily oriented line.

The above pictorial treatment works nicely for the intersection of a horizontal line with a curve. Now consider the intersection of an arbitrarily oriented line with a curve, as illustrated in figure (\ref{fig:newtonsIntersectionAnyOrientation}). Here it is useful to setup the problem algebraically from the begining. Our problem is really still just that of finding the intersection of two lines. The curve itself can be considered the set of end points of the vector

\begin{aligned}\mathbf{r}(x) = x \mathbf{e}_1 + f(x) \mathbf{e}_2,\end{aligned} \hspace{\stretch{1}}(3.3)

for which the tangent direction vector is

\begin{aligned}\mathbf{t}(x) = \frac{d\mathbf{r}}{dx} = \mathbf{e}_1 + f'(x) \mathbf{e}_2.\end{aligned} \hspace{\stretch{1}}(3.4)

The set of points on this tangent, taken at the point x_0, can also be written as a vector, namely

\begin{aligned}(x_0, f(x)) + \alpha \mathbf{t}(x_0).\end{aligned} \hspace{\stretch{1}}(3.5)

For the line to intersect this, suppose we have one point on the line \mathbf{p}_0, and a direction vector for that line \hat{\mathbf{u}}. The points on this line are therefore all the endpoints of

\begin{aligned}\mathbf{p}_0 + \beta \hat{\mathbf{u}}.\end{aligned} \hspace{\stretch{1}}(3.6)

Provided that the tangent and the line of intersection do in fact intersect then our problem becomes finding \alpha or \beta after equating 3.5 and 3.6. This is the solution of

\begin{aligned}(x_0, f(x_0)) + \alpha \mathbf{t}(x_0) = \mathbf{p}_0 + \beta \hat{\mathbf{u}}.\end{aligned} \hspace{\stretch{1}}(3.7)

Since we don’t care which of \alpha or \beta we solve for, setting this up as a matrix equation in two variables isn’t the best approach. Instead we wedge both sides with \mathbf{t}(x_0) (or \hat{\mathbf{u}}), essentially using Cramer’s method. This gives

\begin{aligned}\left((x_0, f(x_0)) -\mathbf{p}_0 \right) \wedge \mathbf{t}(x_0) = \beta \hat{\mathbf{u}} \wedge \mathbf{t}(x_0).\end{aligned} \hspace{\stretch{1}}(3.8)

If the lines are not parallel, then both sides are scalar multiples of \mathbf{e}_1 \wedge \mathbf{e}_2, and dividing out one gets

\begin{aligned}\beta = \frac{\left((x_0, f(x_0)) -\mathbf{p}_0 \right) \wedge \mathbf{t}(x_0)}{\hat{\mathbf{u}} \wedge \mathbf{t}(x_0)}.\end{aligned} \hspace{\stretch{1}}(3.9)

Writing out \mathbf{t}(x_0) = \mathbf{e}_1 + f'(x_0) \mathbf{e}_2, explicitly, this is

\begin{aligned}\beta = \frac{\left((x_0, f(x_0)) -\mathbf{p}_0 \right) \wedge \left(\mathbf{e}_1 + f'(x_0) \mathbf{e}_2\right)}{\hat{\mathbf{u}} \wedge \left(\mathbf{e}_1 + f'(x_0) \mathbf{e}_2\right)}.\end{aligned} \hspace{\stretch{1}}(3.10)

Further, dividing out the common \mathbf{e}_1 \wedge \mathbf{e}_2 bivector, we have a ratio of determinants

\begin{aligned}\beta = \frac{\begin{vmatrix}x_0 -\mathbf{p}_0 \cdot \mathbf{e}_1 & f(x_0) - \mathbf{p}_0 \cdot \mathbf{e}_2 \\ 1 & f'(x_0) \\ \end{vmatrix}}{\begin{vmatrix}\hat{\mathbf{u}} \cdot \mathbf{e}_1 & \hat{\mathbf{u}} \cdot \mathbf{e}_2 \\ 1 & f'(x_0) \\ \end{vmatrix}}.\end{aligned} \hspace{\stretch{1}}(3.11)

The final step in the solution is noting that the point of intersection is just

\begin{aligned}\mathbf{p}_0 + \beta \hat{\mathbf{u}},\end{aligned} \hspace{\stretch{1}}(3.12)

and in particular, the x coordinate of this is the desired result of one step of iteration

\begin{aligned}x_1 = \mathbf{p}_0 \cdot \mathbf{e}_1 + (\hat{\mathbf{u}} \cdot \mathbf{e}_1)\frac{\begin{vmatrix}x_0 -\mathbf{p}_0 \cdot \mathbf{e}_1 & f(x_0) - \mathbf{p}_0 \cdot \mathbf{e}_2 \\ 1 & f'(x_0) \\ \end{vmatrix}}{\begin{vmatrix}\hat{\mathbf{u}} \cdot \mathbf{e}_1 & \hat{\mathbf{u}} \cdot \mathbf{e}_2 \\ 1 & f'(x_0) \\ \end{vmatrix}}.\end{aligned} \hspace{\stretch{1}}(3.13)

This looks a whole lot different than the original x_1 for the horizontal from back at 2.2, but substitution of \hat{\mathbf{u}} = \mathbf{e}_1, and \mathbf{p}_0 = b \mathbf{e}_2, shows that these are identical.

Intersection of two curves.

Can we generalize this any further? It seems reasonable that we would be able to use this Newton’s method technique of following the tangent to refine an approximation for the intersection point of two general curves. This is not expected to be much harder, and the geometric idea is illustrated in figure (\ref{fig:newtonsIntersectionTwoCurves})

Refining an approximation for the intersection of two curves in a plane.

The task at hand is to setup this problem algebraically. Suppose the two curves s(x), and r(x) are parameterized as vectors

\begin{aligned}\mathbf{s}(x) &= x \mathbf{e}_1 + s(x) \mathbf{e}_2 \\ \mathbf{r}(x) &= x \mathbf{e}_1 + r(x) \mathbf{e}_2.\end{aligned} \hspace{\stretch{1}}(4.14)

Tangent direction vectors at the point x_0 are then

\begin{aligned}\mathbf{s}'(x_0) &= \mathbf{e}_1 + s'(x_0) \mathbf{e}_2 \\ \mathbf{r}'(x_0) &= \mathbf{e}_1 + r'(x_0) \mathbf{e}_2.\end{aligned} \hspace{\stretch{1}}(4.16)

The intersection of interest is therefore the solution of

\begin{aligned}(x_0, s(x_0)) + \alpha \mathbf{s}' = (x_0, r(x_0)) + \beta \mathbf{r}'.\end{aligned} \hspace{\stretch{1}}(4.18)

Wedging with one of tangent vectors \mathbf{s}' or \mathbf{r}' provides our solution. Solving for \alpha this is

\begin{aligned}\alpha = \frac{(0, r(x_0) - s(x_0)) \wedge \mathbf{r}'}{\mathbf{s}' \wedge \mathbf{r}'} = \frac{\begin{vmatrix}0 & r(x_0) - s(x_0) \\ \mathbf{r}' \cdot \mathbf{e}_1 & \mathbf{r}' \cdot \mathbf{e}_2 \end{vmatrix}}{\begin{vmatrix}\mathbf{s}' \cdot \mathbf{e}_1 & \mathbf{s}' \cdot \mathbf{e}_2 \\ \mathbf{r}' \cdot \mathbf{e}_1 & \mathbf{r}' \cdot \mathbf{e}_2 \end{vmatrix}}= -\frac{r(x_0) - s(x_0)}{r'(x_0) - s'(x_0) }.\end{aligned} \hspace{\stretch{1}}(4.19)

To finish things off, we just have to calculate the new x coordinate on the line for this value of \alpha, which gives us

\begin{aligned}x_1 = x_0 -\frac{r(x_0) - s(x_0)}{r'(x_0) - s'(x_0) }.\end{aligned} \hspace{\stretch{1}}(4.20)

It is ironic that generalizing things to two curves leads to a tidier result than the more specific line and curve result from 3.13. With a substitution of r(x) = f(x), and s(x) = b, we once again recover the result 2.2, for the horizontal line intersecting a curve.


Having completed the play that I set out to do, the next logical step would be to try the min/max problem that leads to the Hessian. That can be for another day.

Posted in Math and Physics Learning. | Tagged: , | Leave a Comment »

Notes on Goldstein’s Routh’s procedure (continued again.)

Posted by peeterjoot on March 6, 2010

[Click here for a PDF of this post with nicer formatting]Note that this PDF file is formatted in a wide-for-screen layout that is probably not good for printing.

This continues the Routhian procedure notes from last post.

Polar form example.

The troubles appear to come from when there is a velocity coupling in the Kinetic energy term. Let’s try one more example with a simpler velocity coupling, using polar form coordinates in the plane, and a radial potential. Our Lagrangian, and conjugate momenta, and Hamiltonian, respectively are

\begin{aligned}\mathcal{L} &= \frac{1}{{2}} m \left(\dot{r}^2 + r^2 \dot{\theta}^2 \right) - V(r) \\ p_r &= m \dot{r} \\ p_\theta &= m r^2 \dot{\theta} \\ H &= \frac{1}{{2 m}} \left((p_r)^2 + \frac{1}{{r^2}} (p_\theta)^2 \right) + V(r).\end{aligned} \hspace{\stretch{1}}(4.26)

Evaluation of the Euler-Lagrange equations gives us the equations of motion

\begin{aligned}\frac{d}{dt}\left( m \dot{r} \right) &= m r \dot{\theta}^2 - V'(r) \\ \frac{d}{dt}\left( m r^2 \dot{\theta} \right) &= 0.\end{aligned} \hspace{\stretch{1}}(4.30)

Evaluation of the Hamilitonian equations \partial_p H = \dot{q}, \partial_q H = -\dot{p} should give the same results. First for r this gives

\begin{aligned}\frac{1}{{m}} p_r &= \dot{r} \\ -\frac{1}{{m r^3}} (p_\theta)^2 + V'(r) &= -\dot{p}_r.\end{aligned} \hspace{\stretch{1}}(4.32)

The first just defines the canonical momentum (in this case the linear momentum for the radial aspect of the motion), and the second after some rearrangement is

\begin{aligned}m r (\dot{\theta})^2 - V'(r) &= \frac{d}{dt}\left( m \dot{r} \right),\end{aligned} \hspace{\stretch{1}}(4.34)

which is consistent with the Lagrangian approach. For the \theta evaluation of the Hamiltonian equations we get

\begin{aligned}\frac{p_\theta}{m r^2} &= \dot{\theta} \\ 0 &= -\dot{p}_\theta.\end{aligned} \hspace{\stretch{1}}(4.35)

The first again, is implicitly, the definition of our canonical momentum (angular momentum in this case), while the second is the conservation condition on the angular momentum that we expect associated with this ignorable coordinate. So far so good. Everything is as it should be, and there’s nothing new here. Just Lagrangian and Hamiltonian mechanics as usual. But we have two independently calculated results that are the same and the Routhian procedure should generate the same results.

Now, on to the Routhian. There we have a Hamiltionian like sum of p \dot{q} terms over all cyclic coordinates, minus the Lagrangian. Here the \theta coordinate is observed to be that cyclic coordinate, so this is

\begin{aligned}R &= p_\theta \dot{\theta} - \mathcal{L} \\ &= m r^2 \dot{\theta}^2 - \frac{1}{{2}} m \left(\dot{r}^2 + r^2 \dot{\theta}^2 \right) + V(r) \\ &= \frac{1}{{2}} m r^2 \dot{\theta}^2 - \frac{1}{{2}} m \dot{r}^2 + V(r).\end{aligned}

Now, this Routhian can be written in a few different ways. In particular for the \dot{\theta} dependent term of the kinetic energy we can write

\begin{aligned}\frac{1}{{2}} m r^2 \dot{\theta}^2 = \frac{1}{{2 m r^2 }} (p_\theta)^2 = \frac{1}{{2 }} \dot{\theta} p_\theta\end{aligned} \hspace{\stretch{1}}(4.37)

Looking at the troubles obtaining the correct equations of motion from the Routhian, it appears likely that this freedom is where things go wrong. In the Cartesian coordinate description, where there was no coupling between the coordinates in the kinetic energy we had no such freedom. Looking back to Goldstein, I see that he writes the Routhian in terms of a set of explicit variables

\begin{aligned}R = R(q_1, \cdots q_n, p_1, \cdots p_s, \dot{q}_{s+1}, \cdots \dot{q}_n, t) = \sum_{i=1}^{s} \dot{q}_i p_i - \mathcal{L}\end{aligned} \hspace{\stretch{1}}(4.38)

where q_1, \cdots q_s were the cyclic coordinates. Additionally, taking the differential he writes

\begin{aligned}dR &=   \sum_{i=1}^{s} \dot{q}_i dp_i - \sum_{i=1}^s \frac{\partial {\mathcal{L}}}{\partial {q_i}} dq_i- \sum_{i=1}^n \frac{\partial {\mathcal{L}}}{\partial {\dot{q}_i}} d\dot{q}_i- \frac{\partial {\mathcal{L}}}{\partial {t}} dt \\ &=\frac{\partial {R}}{\partial {p_i}} dp_i + \frac{\partial {R}}{\partial {q_i}} dq_i + \frac{\partial {R}}{\partial {\dot{q}_i}} d\dot{q}_i+ \frac{\partial {R}}{\partial {t}} dt,\end{aligned} \hspace{\stretch{1}}(4.39)

with sums implied in the second total differential. It was term by term equivalence of these that led to the Routhian equivalent of the Euler-Lagrange equations for the non-cyclic coordinates, from which we should recover the desired equations of motion. Notable here is that we have no \dot{q}_i for any of the cyclic coordinates q_i.

For this planar radial Lagrangian, it appears that we must write the Routhian, specifically as R = R(r, \theta, p_\theta, \dot{r}), so that we have no explicit dependence on the radial conjugate momentum. That is

\begin{aligned}R &= \frac{1}{{2 m r^2}} (p_\theta)^2 - \frac{1}{{2}} m \dot{r}^2 + V(r).\end{aligned} \hspace{\stretch{1}}(4.41)

As a consequence of 4.39 we should recover the equations of motion by evaluating \delta R/\delta r = 0, and doing so for 4.41 we have

\begin{aligned}\frac{\delta R}{\delta r} = V'(r) - \frac{1}{{m r^3}} (p_\theta)^2 - \frac{d}{dt}\left( -m \dot{r} \right) = 0.\end{aligned} \hspace{\stretch{1}}(4.44)

Good. This agrees with our result from the Lagrangian and Hamiltonian formalisms. On the other hand, if we evaluate this variational derivative for

\begin{aligned}R &= \frac{1}{{2}} m r^2 \dot{\theta}^2 - \frac{1}{{2}} m \dot{r}^2 + V(r),\end{aligned} \hspace{\stretch{1}}(4.43)

something that is formally identical, but written in terms of the “wrong” variables, we get a result that is in fact wrong

\begin{aligned}\frac{\delta R}{\delta r} = m r \dot{\theta}^2 + V'(r) - \frac{d}{dt}\left( -m \dot{r} \right) = 0.\end{aligned} \hspace{\stretch{1}}(4.44)

Here the term that comes from the \dot{\theta} dependent part of the Kinetic energy has an incorrect sign. This was precisely the problem observed in the initial attempt to work the spherical pendulum equations of motion starting from the Routhian.

What variables to use to express the equations is a rather subtle difference, but if we do not get that exactly right the results are garbage. Next step here is go back and revisit the spherical polar pendulum and verify that being more careful with the variables used to express R allows the correct answer to be obtained. That exersize is probably for a different day, and probably a paper only job.

Now, I note that Goldstein includes no problems for this Routhian formalism now that I look, and having worked an example successfully and seeing how we can go wrong, it is not quite clear what his point including this was. Perhaps that will become clearer later. I’d guess that some of the value of this formalism could be once one attempts numerical solutions and finds the cyclic coordinates as a result of a linear approximation of the system equations around the nieghbourhood of some phase space point.

Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

Notes on Goldstein’s Routh’s procedure (continued)

Posted by peeterjoot on March 5, 2010

[Click here for a PDF of this post with nicer formatting]Note that this PDF file is formatted in a wide-for-screen layout that is probably not good for printing.

This continues the Routhian procedure notes started previously.

Now, Goldstein defines the Routhian as

\begin{aligned}R = p_i \dot{q}_i - \mathcal{L},\end{aligned} \hspace{\stretch{1}}(2.16)

where the index i is summed only over the cyclic (ignorable) coordinates. For this spherical pendulum example, this is q_i = \phi, and p_i = m r^2 \sin^2 \theta \dot{\phi}, for

\begin{aligned}R = \frac{1}{{2}} m r^2 \left( -\dot{\theta}^2 + \dot{\phi}^2 \sin^2 \theta \right) + m g r ( 1 + \cos\theta ).\end{aligned} \hspace{\stretch{1}}(2.17)

Now, we should also have for the non-cyclic coordinates, just like the Euler-Lagrange equations

\begin{aligned}\frac{\partial {R}}{\partial {\theta}} = \frac{d}{dt} \frac{\partial {R}}{\partial {\dot{\theta}}}.\end{aligned} \hspace{\stretch{1}}(2.18)

Evaluating this we have

\begin{aligned}m r^2 \sin\theta \cos\theta \dot{\phi}^2 - m g r \sin\theta = \frac{d}{dt} \left( -m r^2 \dot{\theta} \right).\end{aligned} \hspace{\stretch{1}}(2.19)

It would be reasonable now to compare this the \theta Euler-Lagrange equations, but evaluating those we get

\begin{aligned}m r^2 \sin\theta \cos\theta \dot{\phi}^2 + m g r \sin\theta = \frac{d}{dt} \left( m r^2 \dot{\theta} \right).\end{aligned} \hspace{\stretch{1}}(2.20)

Bugger. We’ve got a sign difference on the \dot{\phi}^2 term? But I don’t see any error made.

Simpler planar example.

Having found an inconsistency with Routhian formalism and the concrete example of the spherical pendulum which has a cyclic coordinate as desired, let’s step back slightly, and try a simpler example, artificially constructed

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m (\dot{x}^2 + \dot{y}^2 ) - V(x).\end{aligned} \hspace{\stretch{1}}(3.21)

Our Hamiltonian and Routhian functions are

\begin{aligned}H &= \frac{1}{{2}} m (\dot{x}^2 + \dot{y}^2 ) + V(x) \\ R &= \frac{1}{{2}} m (-\dot{x}^2 + \dot{y}^2 ) + V(x) \end{aligned} \hspace{\stretch{1}}(3.22)

For the non-cyclic coordinate we should have

\begin{aligned}\frac{\partial {R}}{\partial {x}} = \frac{d}{dt} \frac{\partial {R}}{\partial {\dot{x}}},\end{aligned} \hspace{\stretch{1}}(3.24)

which is

\begin{aligned}V'(x) = \frac{d}{dt}\left( - m \dot{x} \right)\end{aligned} \hspace{\stretch{1}}(3.25)

Okay, good, that’s what is expected, and exactly what we get from the Euler-Lagrange equations. This looks good, so where did things go wrong in the spherical pendulum evaluation.


[1] H. Goldstein. Classical mechanics. Cambridge: Addison-Wesley Press, Inc, 1st edition, 1951.

Posted in Math and Physics Learning. | Tagged: , | 3 Comments »

Notes on Goldstein’s Routh’s procedure (setup).

Posted by peeterjoot on March 3, 2010

[Click here for a PDF of this post with nicer formatting]Note that this PDF file is formatted in a wide-for-screen layout that is probably not good for printing.


Attempting study of [1] section 7-2 on Routh’s procedure has been giving me some trouble. It wasn’t “sinking in”, indicating a fundamental misunderstanding, or at least a requirement to work some examples. Here I pick a system, the spherical pendulum, which has the required ignorable coordinate, to illustrate the ideas for myself with something less abstract.


The Lagrangian for the pendulum is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m r^2 \left( \dot{\theta}^2 + \dot{\phi}^2 \sin^2 \theta \right) - m g r ( 1 + \cos\theta ),\end{aligned} \hspace{\stretch{1}}(2.1)

and our conjugate momenta are therefore

\begin{aligned}p_\theta &= \frac{\partial {\mathcal{L}}}{\partial {\dot{\theta}}} = m r^2 \dot{\theta} \\ p_\phi &= \frac{\partial {\mathcal{L}}}{\partial {\dot{\phi}}} = m r^2 \sin^2\theta \dot{\phi}.\end{aligned} \hspace{\stretch{1}}(2.2)

That’s enough to now formulate the Hamiltonian H = \dot{\theta} p_\theta + \dot{\phi} p_\phi - \mathcal{L}, which is

\begin{aligned}H = H(\theta, p_\theta, p_phi) = \frac{1}{{2 m r^2 }} (p_\theta)^2 + \frac{1}{{2 m r^2 \sin^2\theta}} (p_\phi)^2 + m g r ( 1 + \cos\theta ).\end{aligned} \hspace{\stretch{1}}(2.4)

We’ve got the ignorable coordinate \phi here, since the Hamiltonian has no explicit dependence on it. In the Hamiltonian formalism the constant of motion associated with this comes as a consequence of evaluating the Hamiltonian equations. For this system, those are

\begin{aligned}\frac{\partial {H}}{\partial {\theta}} &= - \dot{p}_\theta \\ \frac{\partial {H}}{\partial {\phi}} &= - \dot{p}_\phi \\ \frac{\partial {H}}{\partial {p_\theta}} &= \dot{\theta} \\ \frac{\partial {H}}{\partial {p_\phi}} &= \dot{\phi},\end{aligned} \hspace{\stretch{1}}(2.5)

Or, explicitly,

\begin{aligned}- \dot{p}_\theta &= -m g r \sin\theta - \frac{\cos\theta}{2 m r^2 \sin^3 \theta} (p_\phi)^2  \\ - \dot{p}_\phi &= 0 \\ \dot{\theta} &= \frac{1}{{m r^2 }} p_\theta \\ \dot{\phi} &= \frac{1}{{m r^2 \sin^2 \theta}} p_\phi.\end{aligned} \hspace{\stretch{1}}(2.9)

The second of these provides the integration constant, allowing us to write, p_\phi = \alpha. Once this is done, our Hamiltonian example is reduced to one complete set of conjugate coordinates,

\begin{aligned}H(\theta, p_\theta, \alpha) = \frac{1}{{2 m r^2 }} (p_\theta)^2 + \frac{1}{{2 m r^2 \sin^2\theta}} \alpha^2 + m g r ( 1 + \cos\theta ).\end{aligned} \hspace{\stretch{1}}(2.13)

Goldstein notes that the behaviour of the cyclic coordinate follows by integrating

\begin{aligned}\dot{q}_n = \frac{\partial {H}}{\partial {\alpha}}.\end{aligned} \hspace{\stretch{1}}(2.14)

In this example \alpha = p_\theta, so this is really just one of our Hamiltonian equations

\begin{aligned}\dot{\phi} = \frac{\partial {H}}{\partial {p_\phi}}.\end{aligned} \hspace{\stretch{1}}(2.15)

Okay, good. First part of the mission is accomplished. The setup for Routh’s procedure no longer has anything mysterious to it. Next step for another day.


[1] H. Goldstein. Classical mechanics. Cambridge: Addison-Wesley Press, Inc, 1st edition, 1951.

Posted in Math and Physics Learning. | Tagged: , , | 1 Comment »