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, ``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.