Two finite difference algorithms were implemented in this simulation:
the velocity-Verlet method and the fifth-order
Gear predictor-corrector algorithm. The velocity-Verlet algorithm
involves the following steps:[19,1]

! Step 1: calculate acceleration at time t acc = trapForce + interactionForce ! Step 2: calculate position at time t+dt, velocity at t+dt/2 pos = (pos) + (vel * dt) + (0.5 * acc * dt*dt) vel = (vel) + (0.5 * acc * dt) ! Step 3: calculate acceleration at time t+dt acc = trapForce + interactionForce ! Step 4: calculate velocity at time t+dt vel = (vel) + (0.5 * acc * dt)

In the pseudocode above, ` acc`, ` vel`, and ` pos` are arrays which store the **x**-, **y**-,

and **z**- components
of each ion's acceleration, velocity, and position. Note that the
velocity-Verlet algorithm needs only **3** arrays of dimension .

The fifth-order Gear predictor-corrector method, outlined below, is slightly more complex (but worth it, as we shall see!):

! STEP 1: predict position and up to fourth-order derivatives ! at time t+dt via Taylor expansions... pos = pos+(vel*dt)+(d2*dt2F)+(d3*dt3F)+(d4*dt4F)+(d5*dt5F) vel = vel+(d2*dt) +(d3*dt2F)+(d4*dt3F)+(d5*dt4F) d2 = d2 +(d3*dt) +(d4*dt2F)+(d5*dt3F) d3 = d3 +(d4*dt) +(d5*dt2F) d4 = d4 +(d5*dt) ! STEP 2: evaluate forces on ions at time t+dt acc = trapForce + interactionForce ! STEP 3: "correct" step 1's approximations by a scaled term ! which represents the discrepancy between the second derivative ! obtained from Taylor series expansion and the acceleration ! calculated explicitly from force. del_d2=dt2F*(acc-d2) pos=pos+a0*del_d2 vel=vel+a1*del_d2/dt d2 = d2+a2*del_d2/dt2F d3 = d3+a3*del_d2/dt3F d4 = d4+a4*del_d2/dt4F d5 = d5+a5*del_d2/dt5F

In the pseudocode for the Gear algorithm, ` d2`, ` d3`, ` d4`, and
` d5` are arrays which store the **x-**,
**y**-, and **z**- components of the second, third, fourth, and fifth
derivatives of ions' positions. These derivatives are used in a
Taylor series expansion:

and ` deld_2 `, ` dt2F`, ` dt3F`,
` dt4F`, and ` dt5F`
are factors of for scaling derivative terms.
Using this fifth-order Gear predictor-corrector algorithm requires
**8** arrays of dimension , whereas the
velocity-Verlet algorithm required only **3** such arrays.
The increased memory requirement and complexity of the
Gear algorithm might, at first glance, render it less attractive
than the velocity-Verlet algorithm. However, the Gear algorithm
has an enormous advantage over the velocity-Verlet algorithm:
it requires only * one* calculation of the interaction force per time
step, while
the velocity-Verlet algorithm makes * two* calls to that function at
each update.
Recall that for this long-range interaction problem, the calculation
of interaction force is . Thus, the
speed of the velocity-Verlet algorithm suffers for large **N**.

Moreover, the Gear algorithm can achieve a higher degree of energy
conservation than the Verlet algorithm with a longer time step.
According to Allen,[1] ``based on energy conservation, Gear's
algorithm is about one order higher in accuracy than Verlet's.'' Penning
trap simulations on ` TrapApp` illustrate the superior stability of
the Gear algorithm: one can readily construct a given configuration
and time step such that the velocity-Verlet algorithm will gain energy while
the Gear predictor-corrector algorithm remains stable.

Fri May 12 10:36:01 EDT 1995