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
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 dvThe 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.
(2) -- = v , m -- = f(x) .
dt dt
Starting from a given point (x,v) at time t = 0, the solution to Eq. 2 may be expressed by the parametric equations
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.
dx
x = x(t) , v = v(t), x' = -- .
dt
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
dUthe function U(x) is called the potential energy. Likewise, the kinetic energy is
(3) f(x) = - -- ;
dx
1 2and 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:
(4) K = - m v ,
2
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 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 2The trajectory x(t) which extremizes the action is given by the Euler-Lagrange equation
(5) L(x,x',t) = K - U = - m x' - U(x,t) .
2
d dL dLwhere the derivatives taken with respect to x and x' are understood to be partial derivatives. Substituting Eq. 5 into Eq. 6 yields
(6) -- ------ - -- = 0 ,
dt dx' dx
d dUwhich 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 .
(7) -- (m x') + -- = 0 ,
dt dx
An equivalent formulation may be developed by defining the generalized momentum
dLIn terms of p, the hamiltonian function is
(8) p = ------ = m x' .
dx'
1 2where 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
(9) H(x,p) = px' - L(x,x',t) = - m x' + U(x,t)
2
1 2
= -- p + U(x,t) ,
2m
dp dH dx dHand substituting Eq. 9 into Eq. 10 yields
(10) -- = - -- , -- = -- ,
dt dx dt dp
dp dU dx pwhich is identical in content to Eq. 2.
(11) -- = - -- , -- = - ,
dt dx dt m
x[k+1] = x[k] + h v[k+1/2] ,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!
(12)
v[k+3/2] = v[k+1/2] + h a(x[k+1]) ,
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]) ,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
(13) x[k+1] = x[k] + h v[k+1/2] ,
v[k+1] = v[k+1/2] + (h/2) a(x[k+1]) ;
(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.
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.cdepending on which language you are using. To run the compiled program, give the command
% g77 -o leapint leapint.f
% 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.