PHYSICS 141/241

Winter 2006

Lecture 2: Lagrangian and Hamiltonian Dynamics


2.1 Equations of Motion

The Newtonian equations of motion are second-order in time; thus both positions and velocities are needed to specify the dynamical state of a classical system. The simplest systems are those with one degree of freedom.

Consider a single particle of mass m moving along the x axis under the influence of a position-dependent force f(x). Such a system has the equation of motion

               2
d x(t)
(1) m ------= f(x(t)) .
2
dt

2.2 Phase space

Any second-order equation of the form above can be transformed into a pair of first-order equations by introducing a new variable, the velocity v along the x axis. The transformed equations are

            dx            dv
(2) -- = v , m -- = f(x) .
dt dt
The variables x and v may be visualized as the coordinates of a two-dimensional phase space. At each point (x,v) of phase space, Eq. 2 defines a vector, known as the phase space location.

Starting from a given point (x,v) at time t = 0, the solution to Eq. 2 may be expressed by the parametric equations

   
dx
x = x(t) , v = v(t), x' = -- .
dt
These equations describe a phase space curve passing through the point (X,V) at t = 0. The standard theorems for the existence and uniqueness of solutions to ordinary differential equations imply that there is one and only one phase curve passing through each point.

If f(x) vanishes for some particular X, the phase curve `passing through' (X,0) consists of a single point, known as an equilibrium position. Such points may be stable or unstable; if stable, a point initially in the neighborhood remains there forever, if unstable, it departs in a finite amount of time.


2.3 Conservation of energy

For an isolated system with one degree of freedom, it is always possible to find a function U(x) such that

                     dU
(3) f(x) = - -- ;
dx
the function U(x) is called the potential energy. Likewise, the kinetic energy is
                1    2
(4) K = - m v ,
2
and the sum is the total energy E(x,v) = K + U of the system. Using the chain rule, it is easy to show that the total energy is conserved along any solution of Eq. 2:
            E(x(t),v(t)) = constant.
Thus the phase curves of a system are contours of the total energy.

2.4 Lagrange's equations

Classical mechanics provides an algorithm for obtaining the equations of motion in an arbitrary coordinate system. This algorithm involves defining a function called the lagrangiangeneralized coordinates, the generalized velocities, and possibly the time itself. The equations of motion are then obtained from the principle of least action. which depends on the

For simplicity, this formalism will be illustrated using x as the generalized coordinate and x' as the generalized velocity. The lagrangian is then

                                  1     2
(5) L(x,x',t) = K - U = - m x' - U(x,t) .
2
The trajectory x(t) which extremizes the action is given by the Euler-Lagrange equation
            d    dL     dL
(6) -- ------ - -- = 0 ,
dt dx' dx
where the derivatives taken with respect to x and x' are understood to be partial derivatives. Substituting Eq. 5 into Eq. 6 yields
            d           dU
(7) -- (m x') + -- = 0 ,
dt dx
which is just Newton's law of motion, Eq. 1. The significance of the Lagrangian Action Princliple which provides the unified description of classical and quantum mechanics is expalined in the addendum .

2.5 Hamilton's equations

An equivalent formulation may be developed by defining the generalized momentum

                  dL
(8) p = ------ = m x' .
dx'
In terms of p, the hamiltonian function is
                                         1     2
(9) H(x,p) = px' - L(x,x',t) = - m x' + U(x,t)
2
1 2
= -- p + U(x,t) ,
2m
where x' has been eliminated in favor of p on the second line. Note that numerically the hamiltonian is equal to the total energy. The equations of motion are then
            dp     dH      dx   dH
(10) -- = - -- , -- = -- ,
dt dx dt dp
and substituting Eq. 9 into Eq. 10 yields
            dp     dU      dx   p
(11) -- = - -- , -- = - ,
dt dx dt m
which is identical in content to Eq. 2.

2.6 The Leapfrog Integrator

Many numerical integrators can be applied to a system of coupled ordinary differential eqquations (ODEs), but only some respect the symmetric structure of Eq. 10. One which does is the time-centered leapfrog. Let x[k] = x(k h) be the position at the kth timestep, and let v[k+1/2] = v((k+1/2) h) be the velocity a half-step later. Here h is the timestep which must be held constant throughout the calculation. Then one leapfrog step is
            x[k+1]   = x[k]     + h v[k+1/2]  ,
(12)
v[k+3/2] = v[k+1/2] + h a(x[k+1]) ,
where a(x) = - (1/m) dU/dx is the acceleration. Because this system treats x and v on equal footings, its numerical solution, if computed with infinite precision arithmetic, shares some important properties with Eq. 10. For one thing, both are reversible; equivalently, the numerical result is an exact solution for some other hamiltonian function H'(x,p) near H(x,p). Thus the results describe some hamiltonian, even if not exactly the one specified!

One slight problem with the leapfrog is the need to offset the position and velocity variables by half a timestep. A solution is to split the velocity step:

            v[k+1/2] = v[k]     + (h/2) a(x[k])   ,

(13) x[k+1] = x[k] + h v[k+1/2] ,

v[k+1] = v[k+1/2] + (h/2) a(x[k+1]) ;
this is precisely equivalent to Eq. 12, inasmuch as the relationship between (x[k], v[k+1/2]) and (x[k+1], v[k+3/2]) is the same for both. But when used as a mapping from time k Delta t to time (k+1) Delta t, Eq. 13 is equivalent to starting Eq. 12 with the linear approximation
(14)        v[1/2] = v[0] + (h/2) a(x[0]) + O(h2),
and making a similar linear approximation to extract velocities at integer time-steps. In effect, (a) the solution obtained `jump-starts' from a phase-space point offset in velocity by O(h2) from the specified v[0], and (b) similar errors are made in extracting v[k] at later times. Thus an approximation to the specified hamiltonian is integrated from an approximation of the given initial conditions, and even so the results obtained only approximate the correct trajectory -- and then there are errors due to finite-precision arithmetic.


Laboratory Assignment - First Set

Due date: 1/27/06

For a warmup exercise to the problems you will need to download and modify a simple numerical program. The source code is available in either C (leapint.c) or Fortran 77 (leapint.f). To compile the program, give either of the commands

        % gcc -o leapint leapint.c
% g77 -o leapint leapint.f
depending on which language you are using. To run the compiled program, give the command
        % leapint

The output consists of four columns, listing (1) time, (2) point index, (3) position x, and (4) velocity v. Running the program as supplied yields the trajectory of a point starting at (X,V) = (1,0) in the phase flow defined by the `linear pendulum' or harmonic oscillator, a(x) = -x. The total number of steps taken, number of steps between outputs, and timestep used are determined by the parameters mstep, nout, and dt, respectively; these parameters are set in the main procedure of the program.

Problem 1 

(a)  Modify the statements which set up initial conditions in the main program to produce trajectories starting from the points (2,0) and (3,0). On the (x,v) plane, plot these trajectories (gnuplot or IBM DX) together with the one starting from (1,0).

(b)  Replace the linear pendulum in the accel routine with the `inverse linear pendulum', a(x) = x, and again plot trajectories starting from the points (1,0), (2,0), & (3,0)

(c)  Plot trajectories starting from the same three points for the `nonlinear pendulum', a(x) = - sin(x).

Problem 2 

(a)  Modify the initial conditions to set up n = 100 points in a circle of radius 0.5 centered on the point (0,1), and illustrate the effect of the phase flow of the inverse linear pendulum, a(x) = x, on this circle by plotting these points at several well-spaced times. 

(b)  Now do the same thing for the nonlinear pendulum, a(x) = - sin(x).

Hint: to increase the time interval between successive outputs, use a larger value for the parameter nout. For example, with the given timestep dt = 1/32, setting nout = 32 will output the system state once every time unit.

Problem 3

(a)  Write leapfrog code for the Earth, Mars, and Jupiter orbiting around the Sun. You can also choose Runge-Kutta integration instead. Put the Sun at the origin of your coordinate system. Do some web based research for setting up the initial conditions to have realistic scales of the orbits. Ignore the interactions of the planets in the time integration and keep the Sun fixed at the origin.

(b)  Compare your numerically integrated orbits to the analytic form of the elliptic orbits.

Problem 4

(a)  For each planet, plot the numerically integrated orbit and the analytic elliptic orbit on the same plot. Each planet will have its own plot.

(b)  Animate the orbits simultaneously as the planets are revolving around the Sun. You can choose between Gnuplot and OpenDX.