The Adaptive Verlet Method.
We discuss the integration of autonomous Hamiltonian systems via dynamical rescaling of the vector field (reparameterization of time). Appropriate rescalings (e.g., based on normalization of the vector field or on minimum particle separation in an N-body problem) do not alter the time-reversal symmetry of the flow, and it is desirable to maintain this symmetry under discretization. For standard form mechanical systems without rescaling, this can be achieved by using the explicit leapfrog--Verlet method; we show that explicit time-reversible integration of the reparameterized equations is also possible if the parameterization depends on positions or velocities only.  For general rescalings, a scalar nonlinear equation must be solved at each step, but only one force evaluation is needed. The new method also conserves the angular momentum for an N-body problem. The use of reversible schemes, together with a step control based on normalization of the vector field (arclength reparameterization), is demonstrated in several numerical experiments, including a double pendulum, the Kepler problem, and a three-body problem.
Introduction
. In this article we consider the numerical integration of autonomous
differential equations in R N
d
dt
such as arise in celestial mechanics [22], (classical) atomic and molecular dynamics
[12], and in many other important theoretical and practical situations. Direct integration
(or simulation) is the principal tool for the study of such problems. Fixed stepsize
numerical integration of nonlinear systems (1.1) leads to difficulties, particularly in
the neighborhood of singularities of the vector field. This article describes a family
of efficient adaptive methods for the numerical solution of initial value problems for
equation (1.1) based on normalization of (or other dynamical scaling of) the vector
field.
Recent research [17], [14], [10] has focused on the development of time-stepping
schemes that preserve the underlying geometrical structure of the flow of the system
(e.g. symplecticness in the case of Hamiltonian systems). Regardless of whether a
scheme faithfully preserves available underlying structure, accurate and stable fixed
stepsize numerical integration often requires excessively small timesteps. This is particularly
true when integrating in the vicinity of fixed points and singularities. In
principle, it is possible to incorporate stepsize variation mechanisms based on local
error estimates, but many of the advantages of the structure-preserving schemes appear
to be lost if the stepsize is varied using traditional approaches [4, 3, 7].
We consider an approach based on introduction of a time reparameterization.
Given an appropriate smooth, scalar-valued function
for some m;M , we integrate the differential equations
d
ds
ds
y Department of Mathematics, the University of Kansas, Lawrence, KS 66045, U.S.A.
(whuang@math.ukans.edu). Supported in part by NSF EPSCoR grant no. OSR-9255223.
z Department of Mathematics, the University of Kansas, Lawrence, KS 66045, U.S.A.
(leimkuhl@math.ukans.edu). Supported by NSF grant no. NSF-9303223.
The behavior of trajectories in extended phase space will be highly sensitive to such
reparameterizations, but in terms of the phase space orbits, the reparameterization is
irrelevant. The flow of (1.2) preserves all integral invariants of (1.1). Once the flow
has been reparameterized, (1.2) can be integrated by an appropriate fixed-stepsize
discretization scheme.
A natural choice of R is given by the heuristic R = kfk, with k:k the Euclidean
2-norm. The choice of such a normalization of the vector field can also be viewed
as associating steps in time directly with steps in arclength measured along phase
space orbits. This is not the only choice. Indeed, variants of this reparameterization
idea have been used in various applications, for example basing R on the separation
between particles in a many-particle system with a singular potential [1, 13, 11]. The
techniques of [3, 21] can be viewed in this framework, where R is constructed based on
an embedded local error estimate. Other heuristics would retain approximate constant
arclength of some projection of the solution. A more accurate (and complicated)
control could incorporate an estimate for the orbit curvature at a point. Some related
recent work [18, 19] has used the reparameterization framework to design schemes
that better conserve various integrals such as the energy.
We are particularly interested in systems based on a separated Hamiltonian
It is well known that for large systems of types (1.3)-(1.4) such as one encounters
in molecular dynamics, the computational cost of simulation is dominated by force
evaluations. These problems are often solved by the explicit leapfrog (Verlet) method
which is second order, symplectic, time-reversible (symmetric), and requires only one
force evaluation per step. For large N-body problems, much recent work has concentrated
on finding methods to handle singularities or widely varying forces without
having to evaluate the full force field at every integration step. By basing these methods
on potential-specific splittings, they can sometimes achieve the effect of a variable
stepsize method (see [20]). This article takes a different approach in rather trying to
adapt the time stepsize directly according to some prescribed rescaling of time.
For (1.3)-(1.4), using corresponds to reparameterization based on the
norm of the gradient of the Hamiltonian, hence according to the density of flow lines
along the solution. In general, reparameterization of (1.3)-(1.4) will destroy the Hamiltonian
structure, but time-reversibility of the flow is easily preserved. The coupling
structure of (1.3)-(1.4) is altered by reparameterization, and the natural generalization
of the leapfrog/Verlet scheme (which retains such features as conservation of angular
momentum and time-reversibility) becomes implicit. For example, the schemes of
[21, 3, 11] all lead to implicit integration strategies requiring the evaluation of forces
\GammarV several times per integration step. The main practical contribution of this article
is to show that this implicitness can be removed by the use of a special discretization
method that treats the reparameterization factor as a new variable to be propagated
together with the physical variables. We call this new method Adaptive Verlet.
Actually, we discuss two variants of the method: one based on differentiation of
the equation defining the reparameterization factor, and the other based on introduction
of an algebraic equation which is then solved in tandem with the equations of
motion by use of an ad hoc discretization scheme. This latter differential-algebraic
approach requires no additional force evaluations at a step and it is our recommended
method for typical applications. We stress that for reparameterizations depending on
and q such as that based on arclength, the additional work of the Adaptive Verlet
method over standard Verlet consists only of the solution of a single scalar quartic
polynomial equation at each timestep. Such an efficient implementation would also
be possible for many other choices of the rescaling; in fact, for the special case where
the rescaling is based on either positions or momenta only, the equations become fully
explicit. Control of minimum and maximum stepsizes is also easily incorporated in
our framework.
The methods are implemented for various simple nonlinear models which are
meant to be representative of typical autonomous systems: a double pendulum, the
Kepler problem, and a three-body problem and are found to be robust and efficient
in our experiments.
2. Background. For smooth f , the initial value problem (1.1), u(t 0
at least a local solution through any given u 0 ; we will assume that these solutions
exist for all time. Associated with (1.1) is thus a flow-map OE[t which takes a
given point u 2 R N to the point u  is the solution of (1.1),
. The (forward) trajectory passing through a given point is the set
1g.
There are two ways in which trajectories are unnatural for describing the solutions
of autonomous differential equations. First, because (1.1) is autonomous, it
follows that OE[t (In the sequel, we will write OE[t] for OE[0; t].)
This means that the inclusion of the initial time t 0 in the description of solutions is
inherently arbitrary. Second, if we rescale the vector field solving instead
d
dt
the solution u is unchanged except for a rescaling of time. For these reasons, it is
generally more natural to replace the trajectory through u 0 by the orbit of u 0 defined
as the set fOE[t](u 0 1g. The orbits are scale and translation invariant.
A one-step method for solving (1.1) is a computable approximation to the flow
map, i.e. it is a map /[\Deltat] of R N such that k/[\Deltat](u 0
For a given stepsize \Deltat, the numerical (forward) orbit of a point u 0 is the discrete set
:g. The method is said to be pth order accurate if
fixed and K ? 0 is independent of \Deltat.
The implicit midpoint method defined for (1.1) by
is a second order accurate scheme. A popular method for integrating the separable
system (1.3)-(1.4) is the leapfrog/Verlet scheme:
\Deltat
\DeltatrV (q n+1
4 W. HUANG AND B. LEIMKUHLER
This scheme is explicit, requires only one force evaluation per step, and is 2nd order
accurate.
For the Hamiltonian system (1.3)-(1.4), time-reversibility of the flow is a consequence
of the fact that T
is an even function. This means that
integrating (1.3)-(1.4) from some starting point \Deltat units forward in time, then changing
the sign of the momentum and integrating another \Deltat units returns us to the
starting position. The flow map is also symplectic, meaning that the two form dq -dp
is an integral of motion [6, 17]. These geometric symmetries (symplecticness, time-reversal
impose strong restrictions on the possible types of fixed points
and on the stability of motion in the neighborhood of fixed points.
Recent research on timestepping schemes has concentrated on finding numerical
methods /[\Deltat] which share the time-reversibility and/or symplecticness of the flow.
For example, both the leapfrog and implicit midpoint methods are time-reversible.
They are also both symplectic [6, 16]. The implicit midpoint method (and its Gauss-Legendre
generalizations) preserves all quadratic integral invariants [5], while the
leapfrog scheme (and the broader family of Lobatto IIIa-b pair Partitioned Runge-Kutta
methods) can be shown to preserve the angular momentum [17]. When a
constant stepsize is used, these methods both show remarkable long term stability
with respect to energy error, except in the presence of large forces.
Symplectic (or time-reversible) methods exhibit linear error growth in periodic or
integrable cases [3]. The composition of any two symplectic (time-reversible) maps
is again symplectic (time-reversible), so it is natural to consider varying the stepsize
from step to step. Stepsize selection strategies add a layer of complexity to the
discrete dynamical system, since the stepsize then becomes a function of the dynamic
variables [7]. Experiments with standard variable stepsize integration heuristics have
been disappointing, for example, the general stability of energy and the linear growth
of errors for integrable systems are lost [16, 3]. Time-reversible variable stepsize
methods used in conjunction with underlying symplectic schemes as considered by
Calvo and Hairer [3] have been observed to recapture some of the favorable features
of symplectic/reversible schemes. In this article, we are mainly concerned with the
treatment of generic (nonintegrable) dynamics, but we also observed the linear growth
of errors of [3] for the Kepler problem with our reversible scheme.
3. Singularities and normalization of the vector field. The interactions
of a pair of bodies under a Coulomb potential r is the distance
between the particles) can be regularized by introduction of a new coordinate system
which removes the singularity. In the N-body problem, the two-body interactions
can be similarly regularized [1, 13, 15] but the singularity representing the collision
of three or more particles cannot be eliminated. Similar problems arise with other
singular potentials. This means that some sort of timestep contraction is needed in
the neighborhood of the singularities in order to resolve the behavior accurately. Even
without singularities, it is unlikely that a single, fixed timestep will be appropriate
along a whole orbit.
An insight into the selection of timestep is obtained by noting that for many
problems the norm of the gradient of the Hamiltonian (and hence the norm of the
vector field) is large at the places where the time scale changes rapidly and the fixed
stepsize integrator performs poorly. The norm of the gradient of H measures the
fastest directed rate of change of the Hamiltonian, hence the fixed stepsize integrators
appear to fail at points where many orbits get "squeezed together." To put it another
way, the fixed stepsize integrator may develop instabilities at points where the density
Fig.5a
Fig.5b
Fig.5d
Fig.5c
Fig. 1. Fixed stepsize Verlet method for problem (3.5) with t
of flowlines of the system is large.
As an example of a simple singularity system, we consider the "bond problem"
with Hamiltonian
The orbits of this system in phase space consist of two lobes of closed curves separated
by the line sample orbit is shown in Fig. 1a. Although this orbit is quite
smooth, the timescale changes dramatically in the neighborhood of the origin and the
path from the lower left extreme point on the orbit to the upper left point is traversed
extremely rapidly. In Fig. 1b, this rapid variation appears as an exceptionally steep
derivative of p.
We applied the fixed-stepsize Verlet scheme to integrate the bond problem and
the results for several orbits are shown in Figs. 1c-d with a single fixed stepsize
0:0001. The integration proceeds without difficulty, even for very high energy
levels, as long as the particle stays roughly to the right of the vertical line passing
through the fixed point . Nearer the origin, the curves break-up at relatively
low energy levels and the integrator becomes unstable (Fig. 1d).
Another type of difficulty occurs when the equations of motion possess isolated
near-singular points such as in the double pendulum with very unequal masses. In this
case, cusp-like turning points appear in the solution which can create difficulties for
fixed stepsize integrators. Although a diagonal rescaling (weighting) of the variables
can sometimes remove the near-singular behavior, this is not generically feasible for
complicated problems, so some mechanism is needed to improve the situation.
Now consider a local reparameterization of time according to (1.2) with
kf(u)k. This normalizes the vector field and can be viewed as a control of the arclength
measured along phase space orbits traversed per unit time. A similar heuristic
can be viewed as normalizing the arclength traversed along a trajectory
6 W. HUANG AND B. LEIMKUHLER
per unit time.
Although this article focuses on vector field normalization, the method could be
applied to work with many other choices of R. For example, for particle motion in
central forces, we may typically assume that the closest two-body approach determines
the largest force in the problem. We can introduce the control
min , where
r min represents the smallest interparticle separation and ff is a positive constant. This
idea is suggested in [1, 13, 15, 11].
For the normalized vector field, it is reasonable to use a fixed stepsize integration,
however, it is sometimes necessary to limit the variation of R to avoid excessively small
steps in the neighborhood of singular points. It is easy to incorporate a maximum
and minimum stepsize. For example, for
Given a minimum time stepsize \Deltat min and a maximum time stepsize \Deltat max , and
fictive timestep \Deltas, we set \Deltas=\Deltat min , then we have that
so \Deltat will be bounded between about \Deltat min and \Deltat max .
One drawback of the normalization which becomes evident in the case of separated
Hamiltonians (1.3)-(1.4) is that the normalized vector field is not separable. This
means that explicit schemes, particularly leapfrog/Verlet, cannot be directly applied.
In the next section, we will show how the Verlet method can be modified to provide
an efficient semi-explicit approach to variable stepsize integration. Below, we first
show the effectiveness of combining reparameterization and time-reversible integration
procedures with an example.
Example 3.1: the double pendulum We consider the double pendulum problem
with Hamiltonian
is the gravity constant, m 1 and m 2 are the masses, and the
positions of the bobs are related to ' 1 and ' 2 by
The equations of motion are
To make the problem challenging, we gave the two bobs very different masses (m 1 !!
Because the Hamiltonian is nonseparable, we used the symplectic, time-reversible
implicit midpoint method. A fixed step size method has little difficulty in resolving
x
y
Fig. 2b
x
y
Fig. 2a
timestep
|H-H0|
Fig. 2c
timestep
|H-H0|
Fig. 2d
Fig. 2. The implicit midpoint scheme with fixed stepsize for the double pendulum problem.
are used in the computation. (a),(c)
the solution in the smooth regions, but may fail near the cusps, with the result that
error can increase dramatically in jumps. This behavior is illustrated in Fig. 2. For
fixed stepsize of about 10 \Gamma4 succeeds to resolve the motion
near cusps. When a stepsize of 10 \Gamma3 is used, however, the results are instability in
the energy and poor resolution (especially on long time intervals). The remarkable
stability of energy observed for the smaller value of \Deltat is a hallmark of symplectic (also
time-reversible) methods and is partly explained by the "backward error analysis" of
symplectic methods (see [10], [17]). The cusp points occur where
At these points, in the case of m 1 !! m 2 , the vector field is nearly singular, so we
conclude that the norm of the vector field is large in their vicinity.
For were unable to integrate the double pendulum with stepsizes
larger . The fixed stepsize method is extremely wasteful in terms of
computation, since very large stepsizes are feasible in the smooth regions.
The results of integrating the reparameterized equations for on the
time interval [0; 50] are summarized in Fig. 3. The scheme places many more points
in the neighborhood of the turning points, but uses large timesteps on the smooth
segments. We incorporated a stepsize minimum of \Deltat to restrict the smallest
steps taken in the vicinity of the cusp points. About 380; 000 steps of integration
were needed for this problem which would have required approximately 50; 000; 000
fixed steps for similar accuracy.
4. Semi-explicit integration and the Adaptive Verlet method. In this
section, we turn our attention to the mechanical system with separated Hamiltonian.
8 W. HUANG AND B. LEIMKUHLER
time
|H-H0|
Fig. 4c
time
timestep
Fig. 4d
x
y
Fig. 4b
x
y
Fig. 4a
Fig. 3. The implicit midpoint scheme with the normalized vector field for the double pendulum
are used in the
computation.
The modified equations of motion (1.3)-(1.4) reparameterized in time are
dq
ds
=ae
dp
ds
p). The system is no longer solvable by the standard Verlet method.
The second order Lobatto IIIa-b partitioned Runge-Kutta formula [8] can be
viewed as a generalization of the Verlet method. Applied to (4.7)-(4.8), the equations
are (with
\Deltas
For starting (and to view the scheme as a map of phase space), one may also replace
(4.10) by the two steps
\Deltas
\Deltas
This discretization has the following important features [8]: (i) it reduces to Verlet
it is time-reversible, (iii) it is symplectic if the problem is Hamiltonian
(i.e. if preserves the angular momentum. (Higher order
Lobatto pairs retain these features, but we restrict ourselves to the second order
method in this article.)
Equations (4.9)-(4.10) are implicit, requiring the solution of a nonlinear system
at each step, but the fact that R is scalar can be used to develop a relatively efficient
Newton solver based on rank-one updates. However, still more efficient schemes are
possible.
4.1. Differentiating the control. In this subsection, we describe an approach
based on introducing the stepsize control as a new variable with its own differential
equation. This leads, after discretization, to an algorithm requiring the evaluation of
the Hessian of the potential energy function-or, in a finite difference approximation,
two force evaluations per step. In the next subsection we will see that this scheme
can be further improved.
Differentiating the equation ae = R(q; p) with respect to the reparameterized time
s yields
d
ds
ds
ds
\Theta (r q R) t M
The idea is then to treat the coupled system in (q; p; ae). We will illustrate for the
choice
similar scheme is obtained if upper and lower bounds are built-in to R as described
above.)
R M \Gamma2 p and r q
d
ds
A discretization for the resulting coupled system is obtained by partitioning the
variables as q and (p; ae) and using the second order Lobatto IIIa-b pair:
ae n+ 1p n+ 1(4.11)
ae n+ 1= ae n\Gamma 1+ \Deltas
\Deltas
\Deltas
ae n+ 1(4.14)
(4.12) can be rewritten as
ae n+ 1= a \Deltas
where
Then we get
(ae n+
or, after multiplying by ae 2
(ae
\Deltas
The observation from (4.13) that
ae n+ 1p n+ 1= ae n+ 1p n\Gamma 1\Gamma \DeltasrV (q n )
ae n+ 1ae n\Gamma 1!
is a linear function of ae n+ 1allows us to conclude that (4.15) is a cubic polynomial in
ae n+ 1For starting, we may use the formulas
ae n+ 1= ae n
ae n+ 1;
are given and
The convergence and time-reversibility of this scheme follow directly from the
corresponding features of the underlying Lobatto scheme. The introduction of the
new scalar variables and the fact that the scheme is no longer symplectic do not
destroy the preservation of angular momentum, since, assuming only central forces,
n+1 \Theta p i
n+1 \Theta (p i
n+1 \Theta p i
ae n+ 1p i
n \Theta (p i
In our experiments, we did not use precisely the method described here, but
instead replaced (4.12) by the slightly different midpoint formula
ae n+ 1= ae n\Gamma 1+ \Deltasg(q
The resulting scheme has the same convergence order (two) and the same (linearized)
stability as (4.11)-(4.14), and we observed no evident degradation of performance for
the cases we looked at.
In many cases, the computation of the Hessian matrix V 00 (q) will be an intolerable
expense. In these situations, we suggest using the approximation
for some appropriately chosen ffl. This modification changes only the definition of the
time reparameterization and not the phase space dynamics. It requires two computations
of the forces (at each time step). It may be useful for efficient implementation
to observe that these forces are computed at two very close points in configuration
space; this means that if the forces are evaluated based on lists of nearby interacting
particles, the same list can be used for both calculations.
4.2. A coupled differential-algebraic system. Our preferred approach is to
view the reparameterized equations
dq
ds
=ae
dp
ds
ae
as a coupled differential-algebraic (DAE) system [2]. The method of the previous
section can be viewed as a regularization based on differentiation of the constraints,
and there is a potential for drift effects to degrade the numerical results in long time
interval integration. Since long time interval simulations are precisely the applications
where time-reversible methods can be expected to make a substantial difference, it
may be important to correct this situation. Moreover, we would like to find a method
that does not require computation of the Hessian matrix or more than one force
evaluation per step.
For these reasons, we propose instead the following direct discretization of the
coupled DAEs:
ae
\DeltasrV
ae n+ 1+ ae n\Gamma 1= R(qn
This method requires only a single evaluation of \GammarV at each step.
As for the scheme derived above based on differentiating the control, we can
reduce the implicitness to a single scalar equation for ae n+ 1(ae n+ 1+ ae n\Gamma
for normalized dynamics,
it is easily seen that, after multiplying both sides of (4.23) by ae 2
n+ 1, we are left
with a quartic polynomial in ae n+ 1. In the more general situation of an arbitrary
parameterization R, we would still obtain only a scalar nonlinear equation to solve
for ae.
For the particular situation the computation of ae n+ 1becomes completely
explicit:
ae n+ 1= 2R(qn depends only on p, one can simply reverse the roles of p and q in the staggered
discretization to develop an explicit method.)
The discretization (4.19)-(4.21) is an index-1 differential-algebraic integrator [2],
but not of a standard form. It is easily seen to be second order accurate. We analyzed
the linearized stability of the scheme to confirm the boundedness and convergence of
the iteration under standard assumptions (smoothness of the vector field and boundedness
of the exact solution). In particular, a condition for the stability analysis to go
through is that ae remains strictly bounded away from zero along the exact solution
of interest.
By an argument precisely analogous to that given in x4.1, it is easily seen that
the differential-algebraic method conserves the angular momentum. The method is
clearly time-symmetric.
Now set \Deltat
ae n+ 1. Solving (4.19) for p n+ 1and substituting in (4.20), we
obtain:
\Deltat n
\Deltat n
\Deltat
\Deltat
Thus we can view the method as a traditional variable stepsize version of the St-ormer
rule.
For starting the differential-algebraic scheme from given values (q
could employ the differentiated constraint approach of the previous section. We have
instead used the following starting step with good success:
\Deltas
ae 1= R(q
5. Numerical Experiments. In this section we give several numerical experiments
using the new schemes. In each case, we integrated the normalized vector field
without any minimum stepsize control; a maximum stepsize of \Deltat used. All
experiments were conducted in IEEE floating point arithmetic using double precision
(machine Excellent results were also obtained for the "bond"
problem (3.5).
The use of normalization in conjunction with integrating large N-body problems
is not typically appropriate, since effects of separated strong contributions to the force
field are unreasonably compounded in the determination of the stepsize. Outside of
the well-studied gravitational or purely electrostatic N-body problem with pairwise
potential 1=r (see [1, 13]), we are not aware of any systematic study of appropriate
time transformations for numerical integration. For the small problems considered in
our experiments, normalization of the vector field appears to be a suitable choice.
5.1. The Kepler Problem. The Kepler problem with Hamiltonian
long been a standard benchmark for numerical integration. For low eccentricities,
the orbits can be integrated without difficulty using the standard fixed-stepsize Verlet
method and one recovers with good fidelity all relevant features of the flow. A notable
feature is that, in contrast to standard results for numerical integration, error growth
using a symplectic method such as Verlet is only linear in time.
Our goals for the variable stepsize method are to obtain a method that is robust
with respect to variation of the eccentricity e and which uses small stepsizes only
as orbits pass near the origin (thus improving efficiency). We also demonstrate the
similar behavior of three different methods: the method based on differentiation of
the control, the same method with a finite difference approximation to evaluate the
Hessian, and the differential-algebraic system.
In

Figure

4, we compare the behavior of the three adaptive approaches; all give
very similar results in this experiment. In Figure 5, we see that the adaptive method
easily outperforms the fixed stepsize method. Observe from Fig. 5a that the energy
errors (maximum over the time interval) recorded for the adaptive method are essentially
constant, independent of the eccentricity e. (In fact, the relative energy errors
decrease slightly with increasing e.) The number of steps for various energy errors
(shown in Fig. 5b) is several orders of magnitude smaller than for fixed stepsize. Fig.
5b also shows second order of convergence of the adaptive method, the same as for
the fixed stepsize version.
A more difficult case is shown in Fig. 6, where we used the adaptive scheme to
compute a high eccentricity
14 W. HUANG AND B. LEIMKUHLER
|H-H_0|
Fig.2a
|H-H_0|
Fig.2b
|H-H_0|
Fig.2c
time steps
time
stepsize
Fig.2d
time steps
time
stepsize
Fig.2e
time steps
time
stepsize
Fig.2f
Fig. 4. Comparison of the three adaptive methods for the Kepler problem. t f
0:999 and \Deltas = 0:01 are used. Figs. (a) and (d) for the differentiation formula, (b) and (e) for the
difference formula and (c) and (f ) for the algebraic formula.
adaptive
fixed
1-e
Fig.3a
adaptive
fixed
number of steps
Fig.3b
-0.2
Fig.3c (fixed stepsize)
-0.2
Fig.3d (adaptive)
Fig. 5. Comparison of Adaptive Verlet (differential-algebraic) method and non-adaptive Verlet
methods for the Kepler problem. t f
(plotted every 100 steps). (a) \Deltas = 0:01 and
are used. (b) are used. (d) \Deltas = 0:01 and
are used.
|H-H_0|
Fig.4a
-5
Fig.4b
Fig. 6. Adaptive Verlet (differential-algebraic) method for the Kepler problem. t
are used in the computation. Plotted every 100 steps.
5.2. The three-body problem. As a final example, we looked at a classical
problem that is well known to present a severe challenge for standard numerical
methods: the three-body problem for the gravitational (or electrostatic) interactions
of three point particles. This problem has recently been treated in detail by [15] in
a thorough study of the classical dynamics of Helium, since insight into the quantum
mechanics of the atom can be obtained from an understanding of the bifurcations of
phase space structures such as periodic orbits. There the authors point out, "Quan-
titative classical studies are complicated by the multidimensionality, the nonsepara-
bility, and the unbound character of the motion which is governed by singular and
long-range potentials." The problem is also something of a benchmark in celestial
mechanics.
As a model calculation, we sought to resolve the formation of a binary pair from a
triple following a three-body near-collision. In our experiments, we did not regularize
the two-body collisions as suggested in [15, 1, 13], so somewhat more efficient integration
was probably possible, however, it was the three-body approach that was the
stumbling block in our experiment. In [15], the introduction of a "fictive time"
(time reparameterization) by are the separations of
the two atoms from the frozen nucleus of Helium) is used to improve the behavior of
the dynamical simulation.
Using unit masses and normalizing the gravitational constants, we assigned the
initial positions and momenta
The results of applying the adaptive method with \Deltas = 0:5, \Deltas = 0:1, and
are shown in Figs. 7a,b,c, respectively. The latter two of these integrations
produce nearly identical orbits, and since the timesteps used (Figs. 7d,e,f) are an
order of magnitude smaller in the third case than in the second, we conclude that
with a high probability we have successfully computed the true motion.
Notice that there is a qualitative change in the behavior of the timestep after
about t=3.3. At about this time, the particles are all located in a small region about
y
Fig. 7a
y
Fig. 7b
y
Fig. 7c
time
timestep
Fig. 7d
time
timestep
Fig. 7e
time
timestep
Fig. 7f
Fig. 7. The Adaptive Verlet method on the three-body problem (a),(d)
0:1, and (c),(f) \Deltas = 0:01.
the point (0.25, 2.5) and the forces become very large. One particle is then ejected
and the other two form a binary pair. Very small stepsizes are needed to resolve the
close three-body interaction. Much larger stepsizes can be used to resolve the motion
of the binary pair.
The qualitative change is also illustrated in the energy plots for the three calcu-
lations, shown in Fig. 8.
The average stepsizes used in the three variable stepsize simulations were 0.012,
0.00081 and 0.000073.
We attempted to solve the problem using fixed stepsizes
and . The disastrous results are shown in Figs. 9a,b,c. In all cases, the
integration became completely unstable in the vicinity of the three-body approach.
The energy errors, which were extremely small in the early going, rose dramatically
and suddenly with this event. Fig. 10 shows a close-up of the early behavior for the
case
Only for stepsizes smaller than about is the orbit successfully integrated
with the fixed stepsize method, as shown in Fig. 11.
6. Conclusion. It has long been known that the use of time reparameterizations
can improve the behavior of numerical methods for dynamical systems, particularly
when the orbits pass in the neighborhood of singular points of the vector field. In this
article, we have introduced the Adaptive Verlet method, a time-reversible integrator
for reparameterized N-body problems, and shown it to be both robust and efficient
in numerical experiments. One limitation of this approach is that we do not yet have
a precise (quantitative) understanding of the relationship between the free parameter
\Deltas and the accuracy of the orbit. This is the subject of current work. We are also
investigating efficient higher-order adaptive time-reversible methods, and working to
time
energy
error
Fig. 8
Fig. 8. Energy errors for the Adaptive Verlet method. "Step" refers to \Deltas
apply the method to constrained mechanical systems.

Acknowledgement

Skeel read a preliminary draft of the article and provided
useful comments.
-200200x
y
Fig. 9a
-200200x
y
Fig. 9b
x
y
Fig. 9c
time
energy
error
Fig. 9d
time
energy
error
Fig. 9e
time
energy
error
Fig. 9f
Fig. 9. Fixed stepsize integration of the three-body problem: (a),(d)
time
energy
error
Fig. 10b
x
y
Fig. 10a
Fig. 10. Fixed stepsize integration of the three-body problem, in the vicinity of
(.25, 2.
time
energy
error
Fig. 10b
x
y
Fig. 10a
Fig. 11. Fixed stepsize integration of the three-body problem,
20 W. HUANG AND B. LEIMKUHLER



--R

Integration methods for small N-body systems
Numerical Solution of Initial Value Problems in Differential-Algebraic Equations
Accurate long-term integration of dynamical systems
The development of variable-step symplectic integrators
Stability of Runge-Kutta methods for trajectory problems
A canonical integration technique
Does variable stepsize ruin a symplectic integrator?
Symplectic partitioned Runge-Kutta methods
Matrix Computations
Backward error analysis of numerical integrators and symplectic methods
Building a better leapfrog
Numerical integration methods in moleculardynamics
A numerical experimenter's view of the few-body problem
Explicit canonical methods for Hamiltonian systems
Classical mechanics of two-electron atoms

Chapman and Hall

On the stability of symplectic and energy-momentum algorithms for nonlinear Hamiltonian systems with symmetry
Symplectic integration with variable stepsize
Variable steps for reversible integration methods
Some aspects of chaotic behavior in the solar system
--TR

--CTR
B. Cano , A. Durn, A technique to construct symmetric variable-stepsize linear multistep methods for second-order systems, Mathematics of Computation, v.72 n.244, p.1803-1816, October
B. Cano , A. Durn, Analysis of variable-stepsize linear multistep methods with special emphasis on symmetric ones, Mathematics of Computation, v.72 n.244, p.1769-1801, October
Ben Leimkuhler, An efficient multiple time-scale reversible integrator for the gravitational N-body problem, Applied Numerical Mathematics, v.43 n.1-2, p.175-190, October 2002
n-body problems;variable stepsize methods;hamiltonian systems;leapfrog;verlet;time-reversible methods;symplectic methods