6. Second-Order Dynamical Systems

Astronomy 626: Spring 1997

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.

6.1 Equations of Motion

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

              d x
(1)         m --- = f(x) .

6.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 flow.

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 phase 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 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.

6.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

(3)         f(x) = - -- ;
the function U(x) is called the potential energy. Likewise, the kinetic energy is
                1    2
(4)         T = - m v ,
and the sum is the total energy E(x,v) = T + 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.

6.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 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) .
The trajectory x(t) which extremizes the action is given by the Euler-Lagrange equation
            d    dL     dL
(6)         -- ------ - -- = 0 ,
            dt d xdot   dx
where the derivatives taken with respect to x and xdot are understood to be partial derivatives. Substituting Eq. 5 into Eq. 6 yields
            d             dU
(7)         -- (m xdot) + -- = 0 ,
            dt            dx
which is just Newton's law of motion, Eq. 1.

6.5 Hamilton's equations

An equivalent formulation may be developed by defining the generalized momentum

(8)         p = ------ = m xdot .
                d xdot
In terms of p, the hamiltonian function is
                                            1       2
(9)         H(x,p) = p xdot - L(x,xdot,t) = - m xdot  + U(x,t)
                     1   2
                   = -- p  + U(x,t) ,
where xdot 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.

6.6 The Leapfrog Integrator

Many numerical integrators can be applied to a system of coupled 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 Mississippi half-step later (longer than a New York minute). 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]  ,
            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 (e.g. Saha & Tremaine 1994). 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 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 (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 (t+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(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 O(h^2) 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. Nobody's perfect.



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 -lm
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.

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