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 lagrangian 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 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 Lagrangian Action Princliple provides the unified foundation of the important conservation laws in physics which include energy, momentum, and angular momentum.
(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 motions from the variational principle 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 the Newton formulation as given by Eq. 2.
(11) -- = - -- , -- = - ,
dt dx dt m
Many numerical integrators can be applied to a system of coupled ordinary differential equations (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] ,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.
Tutorials on ODE time integration