PY 502, Computational Physics (Fall 2013)Numerical Solutions of Classical Equations of MotionAnders W. Sandvik, Department of Physics, Boston University1 IntroductionClassical equations of motion, i.e., Newton’s laws, govern the dynamics of systems ranging fromvery large, such as solar systems and galaxies, to very small, such as molecules in liquids and gases(in cases where quantum mechanical fluctuations can be neglected, which is often the case). Inbetween these extremes, Newton’s equations of motion apply, litterally, to ”everything that moves”.Exact analytical solutions of the equations of motion exist only for simple systems, of the typesthat are discussed in elementary classical mechanics courses, and therefore numerical integrationmethods are very important in practice. Here we will discuss some commonly used differentialequation solvers an d use them to study the dynamics of mechanical systems, including ones thatexhibit chaotic dynamics. We will limit the discussion to a single moving body, although themethods can be easily generalized to many-body systems as well—dynamics of many-bod y systemswill be d iscussed later in connection with molecular d y namics simulations.While some of the numerical schemes that we will discuss her e are particularly suitable for integrat-ing classical equations of motions, we will also described m ethods, such as the classic Run ge-Kuttaalgorithm, that are more generally applicable to a large class of ordinary differential equ ations. Thediscussion of systems with chaotic dynamics, although here introduced in the context of classicalmechanics, is also of relevance more br oadly in studies of nonlinear dynamics.2 Basic algorithms for equations of motionConsider a single object (h er e regarded as a point particle) with mass m movin g in one dimension.With its time-dependent position denoted x(t), the differential equation governing its dynamics is¨x(t) =1mF [x(t), ˙x (t), t], (1)where F is the total force acting on the particle, and ˙x and ¨x are the first and second time der ivativesof x. We have indicated that the force may depend on x, ˙x and t. These dependencies typicallycome from from a position dependent static potential, a velocity dependent dampin g (friction),and a time-dep endent driving force, bu t there are other natural posibilities as well, e.g., a positiondependent friction.To study the system numer ically, it is convenient to rewrite the second-order differential equ ation(1) as two coupled first-order equations. Giving the velocity its standard symbol v(t), we have˙x(t) = v(t)˙v(t) =1mF [x(t), v(t), t]. (2)10 10 20 30 4050t-1.5-1-0.500.511.5x(t)0 10 20 30 4050t0.50.60.70.8E(t)∆t=0.01∆t=0.001Figure 1: Time dependence of the position x and the energy E = (1/2)kx2+(1/2)mv2of a harmonicoscillator with k = m = 1 (wh ich gives an oscillation perio d 2π) integrated using the Euler methodwith two different time steps; ∆t= 10−2and 10−3.To integrate this set of equations, we discretize the time-axis as t = t0, t1, . . ., with a constant timestep tn+1− tn= ∆t. The initial values x0= x(t0) and v0= v(t0) are used to start the integrationprocess.2.1 Euler algorithmThe simplest way to advance the time from tnto tn+1is to use the first-order approximation;xn+1= xn+ ∆tvn,vn+1= vn+ ∆tan, (3)where an= F (xn, vn, tn)/m is the acceleration. Clearly, the error made in each step of th isalgorithm is O(∆2t). The method, which is called Euler’s forward metho d, is in general not veryuseful in practice. For example, in systems with no damping or driving force, the energy shouldbe conserved. However, with the Euler method the energy typically diverges with time, whereas inmost higher-order methods the energy errors are bounded. Fig. 1 shows r esults obtained with theEuler method f or the energy and th e position of a harmonic oscillator with k = m = 1 (F = −kx).Even for a time step as as small as 10−3, the energy error is as large as ≈ 5% at t = 50 (correspondingto less than 9 oscillations); it grows faster than linearly with t.There are several ways to proceed to der ive more accurate, higher-order integration schemes; wewill here discuss manipulations leading to a few p ractically useful formulas.2.2 Leapfrog and Verlet algorithmsTo simplify the discussion initially, we will here fir st assume that there is no damping, i.e., theforce and th e acceleration are velocity independent. We begin by writing the Taylor expansion of2v(-1/2) x(0) v(1/2) x(1) v(3/2) x(2) v(5/2) x(3)initial valuesFigure 2: Position and velocity grids used in the leapfrog method. The position is calculated atinteger multiples of th e time step, tn= n ∆t, n = 0, 1, . . ., and the velocity is evaluated at the timestn−1/2= (n − 1/2)∆tbetween these points. The integration starts with given n = 0 values.xn+1= x(tn+ ∆t) to second order in ∆t;x(tn+1) = x(tn) + ∆tv(tn) +12∆2ta(xn, tn) + O(∆3t). (4)Noting that v(tn) + (∆t/2)a(xn, tn) = v(tn+1/2) to order ∆2t, we can rewrite this asx(tn+1) = x(tn) + ∆tv(tn+ ∆t/2) + O(∆3t). (5)We thus need a formula that propagates the velocity on a time grid with points tn+1/2= tn+ ∆t/2,i.e., between the integer-labeled time points tn= t0+ n∆tused for the position. If we use the first-order expansion of the velocity, v(tn+1/2) = v(tn−1/2) + ∆ta(tn−1/2) + O(∆2t), the error remainsO(∆3t) in Eq. (5), bu t we have a problem since this requires the acceleration, and hence the positionx (on which the force depends), on the grid points tn+1/2used only for the velocity. However, itappears intuitively clear, by symmetry, that we should actually use the acceleration at tn+1topropagate the velocity from tn+1/2to tn+3/2, i.e.,v(tn+3/2) = v(tn+1/2) + ∆ta(xn+1, tn+1). (6)The scheme combining Eqs. (5) and (6) is often called the leapfrog method;vn+1/2= vn−1/2+ ∆tan,xn+1= xn+ ∆tvn+1/2. (7)The leapfrog grid is illustrated in Fig. 2.We normally have initial conditions in the form of x0and v0. I n order to start the leapfrog methodwe need v−1/2, which we can get, up to an error ∼ ∆2t, usingv−1/2= v0−∆t2a0+ O(∆2t). (8)The ∆2terror may seem to spoil the O(∆3t) error scaling of the leapfrog method. However, the O(∆3t)error is a sin gle-step error; we will see in Sec. 2.3 that the long-time error, i.e., after integrating alarge number of steps, is O(∆2t), and hence the error in (8) is of appropriate order.It turns out that the position x in the leapfrog