Examples 1253.14.1 The harmonic oscillator
The phase space of a harmonic oscillator with frequencyωand massmwas shown
in Fig. 1.3. Fig. 3.5 (left) shows the comparison between the numerical solution of
the equations of motion ̇x=p/m, ̇p=−mω^2 xusing the velocity Verlet algorithm
analytical solutions. Here, we chooseω= 1,m= 1, ∆t= 0.01,x(0) = 0, andp(0) = 1.
The figure shows that, over the few oscillation periods displayed, the numerical tra-
jectory follows the analytical trajectory nearly perfectly. However, if the difference
|xnum(t)−xanalyt(t)|between the numerical and analytical solutions for the position
is plotted over many periods, it can be seen that the solutions eventually diverge
over time (Fig. 3.5 (middle)) (the error cannot grow indefinitely sincethe motion is
bounded). Nevertheless, the numerical trajectory conservesenergy to within approx-
0 2 4
t/T-1
0
1
x(t)32543 32546
t/T(^050000)
t/T
0
1
2
D(
t)(^050000)
t/T
4.999
5.0
5.001
D
E
(t
) [x10-5]
~
Fig. 3.5(Left) Numerical (solid line) and analytical (dashed line) solutions for a harmonic
oscillator of unit mass and frequencyω= 2. The solutions are shown for 0≤t/T≤4, whereT
is the period, and for 32542≤t/T≤32546.(Middle)Deviation ∆(t)≡ |xnum(t)−xanalyt(t)|.
(Right) Energy conservation measure over time as defined by eqn. (3.14.1)
imately 10−^5 as measured by eqn. (3.14.1), and shown in Fig. 3.5(right). It is also
instructive to consider the time step dependence of ∆Edepicted in Fig. 3.6. The fig-
ure shows log(∆E) vs. ∆tand demonstrates that the time step dependence is a line
with slope 2. This result confirms the fact that the global error in a long trajectory is
∆t^2 as expected. Note that if a fourth-order integration scheme hadbeen employed,
then a similar plot would be expected to yield a line of slope 4.
Next, suppose we express the harmonic forceF(x) =−mω^2 xas
F(x) =−λmω^2 x−(1−λ)mω^2 x=Fλ(x) +F 1 −λ(x). (3.14.2)It is clear that ifλis chosen very close to 1, thenFλ(x) will generate motion on a time
scale much faster thanF 1 −λ(x). Thus, we have a simple example of a multiple time-
scale problem to which the RESPA algorithm can be applied withFfast(x) =Fλ(x)
andFslow(x) =F 1 −λ(x). Note that this examples only serves to illustrate the use of the
RESPA method; it isnotrecommended to separate harmonic forces in this manner!
For the choiceλ= 0.9, Fig. 3.6 shows how the energy conservation for fixedδtvaries
as ∆tis increased. For comparison, the pure velocity Verlet result is also shown.
Fig. 3.6 demonstrates that the RESPA method yields a second orderintegrator with