• Continued discussion on simplectic algorithms
• An example of molecular dynamics: Few body planetary motion.

 

Symplectic algorithms

What is symplectic algorithm?
 

In the Verlet algorithm, our basic dynamic variables are the coordinates of the particles. And we have a somewhat inconvenient situation that the initial conditions are specified as the initial positions and initial velocities, while the difference scheme involves position only. As you may have learned from mechanics, it is possible to reformulate the Newtonian mechanics so that both positions and momenta are the basic dynamic variables. This formulation is called Hamiltonian dynamics (which is equivalent to Newtonian dynamics). A Hamiltonian system is completely described by the Hamiltonian function (kinetic energy plus potential energy):
 
 

H(p,q) = K + V
 
together with Hamilton's equations of motion:
 
 
dp_i/dt = - dH/dq_i,    dq_i/dt = dH/dp_i
 
 
where H is considered to be a function of the generalized coordinates qi and generalized momentum p_i. The derivative on H is a partial derivative. The index i labels the degree of freedom, not the particle.
 

Algorithms for solving this set of  Hamilton's equations of motion are symplectic algorithms.
 

Systems for which a symplectic algorithm can be derived:
 

The symplectic algorithm can be derived if the Hamiltonian takes a simple form:
 
 

H = K + V = sum_i p_i²/2m_i + V(q₁, q₂, ...)
 
This is certainly the case if we use Cartesian coordinates and the interactions between particles depend only on the positions of the particles.

With this particular form of Hamiltonian, Hamilton's equations can be written as:
 
 

dp_i/dt = - dV/dq_i = F_i

dq_i/dt = dK/dp_i = p_i/m_i
 
 
This set of equationa reduce to Newton's equations if one eliminates the variable p_i.
 
 
 

Symplectic algorithms:
 
 

Algorithm A
 
 

We use an Euler type approximation for the momentum derivatives, and then a similar Euler algorithm for the position derivatives. However, for the positions, we use immediately the new result for the momenta. Thus, it is not exactly an Euler algorithm. For brevity, we'll use a vector notation:
 
 
p = ( p1, p2, ... )

q = ( q1, q2, ... )

F = ( F1, F2, ... )
 

and assume that all particles have the same mass m. With these notations and assumption, we have:
 
 
p(t + h) = p(t) + hF(q(t))

q(t + h) = q(t) + h/m p(t + h)
 

Clearly, the algorithm is accurate to first order in h.
 
 
Algorithm B
 
 
We do exactly the same thing as in algorithm A, except that we reverse the order of calculation. The new position is found first.
 
 
q(t + h) = q(t) + h/m p(t)

p(t + h) = p(t) + hF(q(t + h))
 

If you compare algorithms A and B, you will find that it is not exactly the same. It is very interesting to note that both algorithm A and B are in fact mathematically equivalent to the Verlet algorithm if you eliminate the momenta from the equations. However, they have better roundoff error control over the plain Verlet algorithm.
 
 
Algorithm C (second-order symplectic algorithm or velocity form of the Verlet algorithm):
 
 
This is the most widely used algorithm that we are really interested.

This algorithm is derived by dividing the step size h into two, for the first half step size h/2, we use the algorithm A, for the second half of the step, we use the algorithm B. In total, we have full step size h. That is, use the formula for algorithm A but set h/2. We get equations relating the quantities at time t to time t + h/2. Then we use the equations for algorithms B with step size again h/2, but the time is for t+h/2 to t + h. We get four equations:
 
 

p(t + h/2) = p(t) + h/2 F(q(t))

q(t + h/2) = q(t) + h/2m p(t + h/2)

q(t + h) = q(t + h/2) + h/2m p(t + h/2)

p(t + h) = p(t + h/2) + h/2 F(q(t + h))
 

Now we eliminate the intermediate quantities evaluated at time t + h/2, we get the final result for the symplectic algorithm:
 
 
q(t + h) = q(t) + h/m p(t) + h2/2m F(q(t))

p(t + h) = p(t) + h/2 [ F(q(t)) + F(q(t + h)) ]
 
 

Note that the coordinates q have to be evaluated first since the new values at time t+h are used to evaluate the forces.

The algorithm C is also second-order in step size h. From the derivation it is not obvious that it is necessarily superior than Verlet algorithm. However, the symplectic algorithm has several advantages:
 

• It is time-reversible
 
• The system can be started naturally, with initial position q0 and initial momentum p0 = mv0
 
• The symplectic algorithms A, B, and C have one important property that share with the original Hamiltonian system---they preserve Poincare invariants.
 
 
This last property of Hamiltonian system is usually discussed in advanced mechanics or dynamical system courses. It is about the phase space (p,q) flow property. A solution with initial condition (p0,q0) of the Hamiltonian system p=p(t,p0,q0), q = q(t,p0, q0) can be viewed as a map in the phase space from the point (p0,q0) to the point (p,q), where the time parameter t is taken as some constant. The map is also called a canonical transformation, or symplectic transformation.

One of the series of Poincare invariants is the phase space volume. Certain region of points of volume V0 can be mapped to a set of new points forming volume V. Invariance means V=V0. This is called Liouville theorem.

Not only is the volume in phase space invariant under the mapping of the dynamics, but also there are other lower dimensional objects (hypersurfaces) which are invariant. The symplectic algorithms preserve this invariant property exactly.

 

 
 

 
 

Two body problem:
 

E.g.: Earth and Moon.

With Earth fixed at the center, the motion of the Moon is exactly solvable.

Three body problem:
 

E.g. Sun, Earth and Moon.

No analytical solution.
 

Numerical solution of the three-body problem:

Consiter a system consisting of the Sun, which is fixed at the center of coordinate system and with mass M, the Earth with mass m and position r, and the Moon with mass m' and position r'.

For simplicity, we assume that both the Earth and the Moon move in the same two-dimensional plane, so that the vectors will be understood as two-dimensional vectors containing just x and y components.

The magnitude of the gravitational force between any pair of objects is:
 

F = G mM/r²

that is, the gravitational force is proportional to the masses of the objects, and is inversely proportional to the distance square. The direction of the force is always attractive and is towards the other object.

The gravitational constant is

G ~ 6.673 x10^-11 m³ kg^-1 s^-2

We have two coupled Newton's equations of motion, one for the Earth and one for the Moon:
 

m d²r/dt²
= - G mMr/r³ - G mm'(r-r')/|r-r'|³

m' d²r'/dt²
= - G m'Mr'/r'³ - G mm'(r-r')/|r-r'|³
 

where r = |r| = sqrt{x2+y2}.

The first term of the both equations on the right-hand side is the attractive force due to the Sun, the second term is mutual interactions between the Earth and the Moon. The vector notation takes care of the magnitude as well as the directions of the forces.
 

 

Computer algorithm development:

Before plugging these equations into a computer, we need to change all the symbols to numbers. If we want to study the actual motions of the Sun-Earth-Moon system, we'd better use the correct constants for G, M, m, and m'. However, we should not do this naively, since they are very large or very small numbers that we might get overflow or underflow in calculations. It is better to normalize the basic units so that the typical quantities are not very far from 1.

In this particular case, for the purpose of exploring the general feature, we can take, say, G=1, M=100, m=1, and m'=0.1. This amounts to take different units than second, meter, and kilometer and to use a particular mass ratios of the objects. With this strategy in mind, it is a straightforward task to implement a numerical integration of the Newton's equations.

For example, if we use the Verlet algorithm, we can follow the motion of the Earth and Moon as:
 

r(t+h)=-r(t-h)+2r(t)
            +h²[-GM r(t)/r³(t)-Gm'(r(t)-r'(t))/|r(t)-r'(t)|³

r'(t+h)=-r'(t-h)+2r'(t)
             +h²[-GMr'(t)/r'³(t)+Gm(r(t)-r'(t))/|r(t)-r'(t)|³
 
 
 

Homework