Page 205 - Numerical Methods for Chemical Engineering

194 4 Initial value problems

grid points, and that the slowest modes have longer wavelengths comparable to the dimen-
sions of the physical domain. Often, the real physical behavior occurs on these long wave-
lengths and the true solution has little variation from one grid point to the next. Thus,
for the true solution, the fast modes with [Re(λ j )] −1 ( t) should have c j ≈ 0, and we
only observe significant activity in these modes, c j = 0, due to the exponential growth of
round-off errors when |µ j | 1. If we use a BDF method with stiff decay, even though we
are not simulating the fast mode dynamics accurately, they do correctly decay to zero in the
true solution. The fast modes therefore do not corrupt the dynamics of interest taking place
on the slow time scales. Therefore, integrators with stiff decay can simulate the dynamics
occurring on the slow time scales of a stiff system using time steps that are very large on
the scale of the fast modes. Essentially, a QSSA is made on the fast mode dynamics.

Symplectic methods for classical mechanics

The integration methods discussed above have been compared in terms of stability and
accuracy; however, no integration method with a finite time step is perfectly accurate. Over
long simulation periods, the discrepancy between the predicted and actual system behavior
may be significant.
This is of importance for applications such as celestial mechanics and molecular dynam-
ics, in which we simulate the motion of a number of interacting particles of masses m α ,
positions r α , velocities v α , with a total potential energy function V (r 1 ,r 2 ,..., r N ). The
motion of each point mass is governed by Newton’s second law of motion

∂
dv α
m α = F α F α =− V (r 1 , r 2 ,..., r N ) (4.183)
dt ∂r α
In the absence of an external potential or dissipation (frictional forces), the total system
energy is constant,

N

1
E = m α (v α · v α ) + V (r 1 , r 2 ,..., r N ) = constant (4.184)
2
α=1
Due to integration errors, this property is generally not satisfied by the numerical trajec-
tory of the system; however, for a special class of integrators, total energy is conserved
(approximately) even over long simulations.
From Noether’s theorem (Arnold, 1989), it is known that the conservation of energy is
related to the invariance of Newton’s equation of motion to time reversal. That is, if we
follow a conservative system for some period of time and then reverse the direction of time,
the system will exactly retrace, in reverse, its previous trajectory. It may be shown that
integration rules that are symplectic (i.e., symmetric with respect to the direction of time)
have favorable energy conservation properties that make them most suitable for simulating
the classical mechanics of conservative systems. A more complete discussion of symplectic
integrators is found in Frenkel & Smit (2001). Here, we merely provide a popular symplectic