AN ENERGY CONSERVING MODIFICATION OF NUMERICAL METHODS FOR THE INTEGRATION OF EQUATIONS OF MOTION

In the integration of the equations of motion of a system of particles, con- ventional numerical methods generate an error in the total energy of the same order as the truncation error. A simple modification of these methods is described, which re- sults in exact conservation of the energy.


INTRODUCTION.
When applied to the motion of a system of particles, conventional numerical methods for the integration of ordinary differential equations only approximately conserve the total energy of the system.The error in the calculated value of the energy is of the same order as the truncation error in the velocities.In previous work [1]- [5], a new class of methods was described, which maximally conserve the constants of motion.
These methods exactly conserve the total energy and linear momentum, and conserve the total angular momentum to at least one higher order than the corresponding conventional methods.
In what follows, our purpose is to show how conventional numerical methods--ex- emplified by the third-order Taylor series and Adams' formulae--can be modified so that exact conservation of energy occurs.This modification simply involves the in- troduction of adjustable, multiplicative parameters, whose values are unity for the conventional case.

EQUATIONS OF MOTION.
The following is a brief description of the equations of motion of a system of n particles, interacting according to a pairwise-additive potential.For more details, see [i] or [5].
Suppose particle i has mass mi, position vector r.

I i I
relates the acceleration a. to the force Fi, given by 1 r. (2.4) where is the potential of interaction.It will be assumed that has the pair- wise-additive form For n particles, equation (2.10) yields a system of second-order ordinary differ- ential equations for the r..1 This system may be used to solve for the and i' at any later time t' t + At, given the r. and v. at time t.
1 1 Conservation of the total energy E occurs because of the existence of the poten- tial .Here, where a-b denotes the scalar product of two vectors a and b.Conservation of energy is expressed by the equation for any two times t and t', with E evaluated along the trajectory.

NUMERICAL METHODS FOR THE INTEGRATION OF EQUATIONS OF MOTION 175 3. CONVENTIONAL NUMERICAL METHODS
A simple example of a conventional approximation method for the numerical solu- tion of equations (2.10) is provided by the truncated Taylor-series formulae where the r'c,i and V'c,i are the calculated values for the r ol and v_ at time t' t + At, and dq d0ji vji 0ji i d0 jirji   When equations (4.1) and (4.2) are used to obtain estimates for r.l and v i, an error AE is made in the total energy, which is given by  (preserving the order of the method) and so that exact conservation of energy occurs.
for the e.. gives, for example, for (4.5), the equation (cf. [i] and [5]) Because the equations of (4.7) are linear except for terms of 0[(At)3], they may be easily solved via the iteration formula &ij At + (V.o+a. -) . ... i (4.9) Iteration to convergence of the e.. guarantees exact conservation of energy in the me tho d.
Higher-order formulae may be obtained directly in the same way as equation (4.7). in equation (4.3), the sum transformed to i <j, and the ij terms set individually to zero.These resulting implicit equations in the e.. are then solved by standard methods, with the first approximations given by equations (4.9).For very high order methods, the extra algebra needed to obtain the e.. is con- siderable, and substantially reduces the relative efficiency of the method.However, it should be noted that conservation of energy guarantees stability in the usual sense (bounded motion), which is always a desirable computational property.

NUMERICAL EXAMPLE.
As an illustration of the affect of the modification described in Section 4, the modified and unmodified forms of the third-order Adams' method are compared numerically on a sample two-dimensional problem involving two particles.
Here n 2, m I m 2 2 (5.1) and the gravitational interaction i (5.2) 12(r12) r12 is used.The initial conditions are chosen so that the center-of-mass of the system is at rest with r12(0) (,0)   (5.3) v12(0) (0,1.63). (5.4) The value of the energy is then E 0.6715500000... (5.5) Because of the form of 12 in (5.2), the exact motion that occurs traces out a closed ellipse with major-axis 2a 1.48909 23855  (5.6)   corresponding to upper and lower bounds on r12 of r> 0.98909 23855 and r< 0.50000 00000.
The motion repeats itself with period equal to 4.0366 15087, The implicit equations of the third-order methods were iterated to a relative con- -8 ergence of I0 A constant step-size of t /80 was used.In order to focus attention on the errors made in the methods, results were obtained at times t which were multiples of the period where the exact solu- tion returns to the initial conditions.Measures of the errors at these points are the error in the calculated value of E, the deviations from zero of dX/dt and , and the deviation of r12 from 1/2.It can be seen that the unmodified Adams' method makes an error in E as well as larger errors in dX/dt and Y, and compares unfavorably with the modified method.Another simple measure of the error for this problem is the number of steps over which a phase error of 180 is made: i.e., the time at which r12 0.985 instead of 0.5.For the unmodified methods, this was about 2800 steps (35T).For the modified methods, at 20000 steps (250T), a phase error of less than 180 had been made.
Programs for the methods are given in the Appendix of [6].

6 )
equations (3.1) and(3.2) is of third-order, since 'r.7 + 0[(At) 4neglect of the succeeding Taylor-series terms.These errors generate an error of 0[(At) 3] in the value of the energy E' calculated using the 'The third-order Adams' method arises via equations (3.1) and (3.2) and the ap- equation (3.8) ' denotes the value of c,ij using the r'c,zj+o.. Equations (3.4), (3.5), and (3.6) also hold when the ij used for the Gij.
Fij obtained from equation (2.9) are 4. ENERGY CONSERVING MODIFICATION OF CONVENTIONAL METHODS.

(4. 7 )
For n particles, equation (4.7) yields a set of implicit, coupled equations in the depend upon the values of the ij' since the bi3 and c,ij ij" For small At the equations of (4.7) are strongly linear in the ... The e ij c,ij' (through only nonlinear dependences on the occur through the bij and the r' ..).In both these cases, the terms involving the e.o occur with coefficients c, 13 lJ proportional to (At)3.(Compare equations (4.1), (4.2), (4.4), and (4.7).)In con- trast the coefficients of the linear terms in ij' namely At (vij + aijAt) ijare of 0[At].
the S.o satisfy equation(4.6).The formulae for the v' Consider the third-order methods of Section 3, with G.*. replacing either the

Table I
gives these quatities for several times t mr.
CThird-order Adams' method modified to give exact energy conservation.

TABLE I .
Comparison of Modified and Unmodified liethods on a Simple Gravitation Problem