(Comp.sys.hp48) Item: 693 by charliep@hpcvra.cv.hp.com. [Charles Patton] Subj: Hamiltonian Method for Geodesics Date: Wed Feb 19 1992 @ Version 1.0 @ _Introduction_ The programs contained in this document comprise a collection of simple-minded utilities for use in computing approximate orbits of Hamiltonian dynamical systems in three space variables. These utilites are the outcome of messing around with ways of computing geodesics on implicitly-defined surfaces embedded in Euclidean space. There is nothing new about the method involved, rather, the nicest thing is how easy it is to "follow your nose" when working in an environment like that of the '48. All suggestions, additions, corrections, insights, commentary, and criticisms are welcome. _NO WARRANTY_ This work is provided on an "as is" basis. Hewlett-Packard Co. provides no warranty whatsoever, either express or implied regarding the work, including warranties with respect to its merchantability or fitness for any purpose whatsoever. _Copyright_ Copyright (C), 1991, Hewlett-Packard Co. Permission to copy all or part of this work is granted provided that the copies are not made or distributed for resale (excepting nominal copying fees) and the _NO WARRANTY_ and this copyright notice are included verbatim. Other permissions can be arranged by contacting the author. _Correspondence_ Correspondence on HamiltonGeodesy should be sent to Charles Patton at one of the following addresses: snail-mail: M.S. 5U-L9 Hewlett-Packard Co. 1000 N.E. Circle Blvd. Corvallis, OR 97330 e-mail: charliep@cv.hp.com (for Internet hosts) hplabs!hpcvrs!charliep (for UUCP hosts) _Overview_ Recall that the Hamilton-Jacobi formulation of classical mechanics starts with the energy function ("Hamiltonian") as a function of position and momentum of the system in question: q,p --> H(q,p). The Hamilton equations then say that the time derivative of the position is the derivative of H with respect to momentum, and the time derivative of the momentum is the negative of the derivative of H with respect to the position: q'=dH/dp ; p'=-dH/dq. The orbits of the Hamiltonian are the integral curves of this system of differential equations, that is, curves q(t),p(t) which satisfy dq/dt=dH/dp and dp/dt=-dH/dq. Perhaps the simplest way to formulate the problem of finding geodesics on a Riemannian space is as a Hamiltonian system. In this case, the Hamiltonian is just the function which assigns to a position and momentum half the squared length (according to the local metric) of the momentum: H(q,p)=0.5 Sum g (q) p p i,j ij i j where g(q) is the metric at q. In this case, the projection onto position space of the orbits of the Hamiltonian are the geodesics. If we are trying to find geodesics on a surface emdedded in Euclidean space, we could find local coordinates for the surface, compute the metric in local coordinates, and then use the above formulation to find geodesics in terms of these local coordinates. This is almost always difficult to do. We will take a different approach here, the key difference being that we are assuming that the surface is given to us as the zero-set of a function F: q --> F(q) defined on the Euclidean space. We can use F to define a potential function on the Euclidean space with a minimum along the desired surface and use the Euclidean metric together to define the Hamiltonian: 2 2 2 H(q,p)=k F(q) +0.5 ||p|| for some "large" constant, k. If we start with any initial q,p with F(q)=0 and p "tangent" to the the surface, the projection onto position space of the resulting orbit will approximate a geodesic on the desired surface. You can think of this as describing a particle attached to the surface by a sliding, frictionless spring. If the spring is very stiff, the particle will follow the path of a geodesic when it is given a shove along the surface. If the spring is not so stiff, the particle will combine near-geodesic motion with oscillations normal to the surface. _Organization and Descriptions of the Programs_ The utilities are arranged in a directory with the two top-level programs first, and the supporting utilities and variables afterwords. The space coordinates are represented by the variables, 'x' 'y' and 'z' (note the lower case.) The corresponding momentum coordinates are 'px' 'py' and 'pz'. SetUpHamiltonian: Given the function kF, on the stack, as an expression in the space variables, 'x' 'y' and 'z', this program creates the Hamiltonian described above and stores it in the variable 'H' for future reference. It then computes the components, 'xdot'...'pzdot', corresponding to the Hamiltonian equations. It then calls on the utility Lie4thOrder to compute a (formal) fourth order polynomial approximation to the solution of the system of Hamilton equations, storing the formulae for the increments as 'x'...'pz'. This takes a loooong time so go out for a coffee break after firing it up. After this finishes, you are ready to plot approximate geodesics on your surface. You can easily customize this program to set up for an arbitrary Hamiltonian in these variables by deleting the first line of the program listed, below, and storing your choice of Hamiltonian (as an expression in 'x'...'pz') in the variable 'H' before running SetUpHamiltonian. N.B.>Insure that 'x' 'y' 'z' 'px' 'py' and 'pz' are formal (no variables N.B.>with those names along the current path.) N.B.> Examine the MODES menu to be sure that the N.B.> the machine is in SYMBOLIC EVALUATION mode ( SYM button.) RunHamiltonian: This program takes initial values for x y z px py and pz on the stack and computes successive points (approximately) on the orbit of the Hamiltonian. It graphs these by lines connecting successive (x,z) pairs (i.e. projection of the orbit on the x-z plane) using the current PPAR values for the horizontal and vertical scaling of the screen. It requires that you have previously run SetUpHamiltonian in the current directory. As part of its initialization, it takes the current value of the stepsize (stored in the variable ) and pre-computes the higher order stepsizes (^2/2, ^3/6, and ^4/24) storing these in 2 3 and 4. As explained below, the default value of is 0.1, but you can use any value you wish as long as you store it before running this program. You should exit this program by pressing the [ENTER] key. When this is done, the program will clean up and leave the last computed values of x...pz on the stack ready for continuation of the orbit, if desired. HamiltonStep: This simple utility takes the values off the stack, stores them in 'x'...'pz', computes the next step and returns the values to the stack. Lie4thOrder: This routine constitutes the heart of this collection of utilities. The idea is that applying the Lie derivative operator d d d d d d L:=x'- + y'- + z'- + px'- + py'- + pz'- dx dy dz dpx dpy dpz to any function, f, corresponds to taking the time derivative of f along the orbits of the Hamiltonian. Thus, the fourth-order Taylor series approximation to the time evolution of f is f+(Lf)t+(LLf)t^2/2+(LLLf)t^3/6+(LLLLf)t^4/24 Using this on each of the coordinate functions, x...pz, thus gives an approximation to the orbits themselves. This program computes this approximation for each of the coordinate functions using as the time-step, t. This is equivalent (in the limit of small time-step) to the fourth-order Runge-Kutte method of approximating solutions to ordinary differential equations, but has the advantage of being much more comprehensible. It also allows you to compute the approximate time evolution of any smooth function, not just the coordinate functions, in exactly the same manner. The resulting approximation is subject to the same kinds of instabilities as the R-K method. In our case, this shows up when the approximation overshoots into a region where the added potential is steep. This causes the next step to overshoot even more in the other direction, quickly leading to huge values for the position and momentum coordinates. For this reason, you should leave at its default value, and if instabilites become apparent, choose a smaller value for the initial velocity. The point here being that while the space-projection of the orbits of the geodesic equation don't depend on the initial velocity, the orbits of our approximation do. They are only asymptotic to the geodesics in the limit of small initial velocity. Moreover, no choice of potential will fix that. You can think of it like a motorcycle driving in a groove. At low velocity, it will follow the groove nicely, but at high velocity, it will "bank" its turns on the wall of the groove, hence the possibility for overshoot becomes more serious. LieDer: This utility takes a function, given as an expression in x...pz, and computes its Lie derivative according to the current values of 'xdot'...'pzdot' which correspond to x'...pz'. The Other Variables: ...4 and 'H' were explained previously. The others are placeholders for values computed during the operations of the main programs.