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 (1) m --- = f(x) . 2 dt

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 dtThe variables

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

x = x(t) , v = v(t) .These equations describe a

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.

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) = - -- ; dxthe function

1 2 (4) T = - m v , 2and the sum is the

E(x(t),v(t)) = constant.Thus the phase curves of a system are contours of the total energy.

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

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

1 2 (5) L(x,xdot,t) = T - U = - m xdot - U(x,t) . 2The trajectory

d dL dL (6) -- ------ - -- = 0 , dt d xdot dxwhere the derivatives taken with respect to

d dU (7) -- (m xdot) + -- = 0 , dt dxwhich is just Newton's law of motion, Eq. 1.

An equivalent formulation may be developed by defining the generalized momentum

dL (8) p = ------ = m xdot . d xdotIn terms of

1 2 (9) H(x,p) = p xdot - L(x,xdot,t) = - m xdot + U(x,t) 2 1 2 = -- p + U(x,t) , 2mwhere

dp dH dx dH (10) -- = - -- , -- = -- , dt dx dt dpand substituting Eq. 9 into Eq. 10 yields

dp dU dx p (11) -- = - -- , -- = - , dt dx dt mwhich is identical in content to Eq. 2.

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

One slight problem with the leapfrog is the need to offset the position and velocity variables by half a timestep. A lazy 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

(14) v[1/2] = v[0] + (h/2) a(x[0]) + O(h^2),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

- Arnold, V.I. 1978,
*Mathematical Methods of Classical Mechanics*. - Binney, J. & Tremaine, S. 1987,
*Galactic Dynamics*, Appendix 1.D. - Saha, P. & Tremaine, S. 1994,
*A. J.***108**, 1962.

Due date: 2/06/97

For the next two 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

% cc -o leapint leapint.c -lm % f77 -o leapint leapint.f -lmdepending 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.

6. (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 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)`.

7. (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.*

Joshua E. Barnes (barnes@galileo.ifa.hawaii.edu)

Last modified: February 1, 1997