Next: Choosing a Time Up: Methods and Theory Previous: System of Units

# Finite Difference Methods

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.

Next: Choosing a Time Up: Methods and Theory Previous: System of Units

Wolfgang Christian
Fri May 12 10:36:01 EDT 1995