Finite Difference Methods

Two finite difference algorithms were implemented in this simulation: the user chooses at run-time between the velocity-Verlet method and the Gear predictor-corrector algorithm. The velocity Verlet algorithm involves the following steps [1][9]:

.
! Step 1: calculate acceleration at time t
      CALL trapForce
          ! returns force due to trap to variable acc
      acc=acc+interactionForce()
! Step 2: calculate position at time t+dt, velocity at t+dt/2
      pos=pos+vel*dt+0.5*acc*dt**2
      vel=vel+0.5*acc*dt

! Step 3: calculate acceleration at time t+dt
      CALL trapForce
      acc=acc+interactionForce()

! Step 4: calculate velocity at time t+dt
      vel=vel+0.5*acc*dt

In the code above, acc, vel, and pos are arrays which store the -, -,
and - components of each ion's acceleration, velocity, and position. Note that the velocity-Verlet algorithm needs only arrays of dimension .

The Gear predictor-corrector method is slightly more complex:

.
! 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
         CALL trapForce
         acc=acc+interactionForce()

! STEP 3: "correct" step 1's approximations via scaled term
! which represents the discrepancy between Taylor-series second
! derivative and the acceleration calculated explicity 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 code for the Gear algorithm, d2, d3, d4, and d5 are arrays which store the , -, and - 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 arrays of dimension , whereas the velocity-Verlet algorithm required only 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, one might expect the speed of the velocity-Verlet algorithm to suffer for large . Indeed, Table 1 shows that run times are comparable for simulations involving a very small number of ions, but the time required by the Verlet algorithm is forty times greater for a simulation of ions.

Once a finite difference algorithm is chosen and implemented, one must determine the optimal time step, , for the simulation. A ``good'' time step is large enough that the system evolves rapidly, yet small enough that the finite difference algorithm is accurate. There is no inherently ``good'' time step- rather, the size of an appropriate time step depends on the system being simulated. In this simulation, the Einstein period is used to scale the time step to a value which is suitable for the simulation. The Einstein period is with the Einstein frequency, given by

where and are the average force and velocity squared acting on each ion with mass .

Allen and Tildesley suggest that the time step, , should be less than the Einstein period by at least an order of magnitude for simulations of liquids [1]. However, Table 2 shows that this ion trapping simulation requires a much smaller time step. A variety of time steps were tested for stability and accuracy. A numerically unstable algorithm tends to ``blow up'' and result in overflow errors; an inaccurate algorithm yields results which are incorrect, but may appear well-behaved. Note that numerical instability results even when is as small as . Setting resulted in stable, accurate results for Penning trap and Paul trap simulations showing both trapping and escaping behaviors.

Using the Einstein frequency to scale the time step was crucial to the program, because it relieves the user of the task of finding an appropriate time step for each simulation and allows the program to use a smaller number of time steps than if an artificially small time step is ``hard-wired'' into the program. Before the time step was implemented as a function of , was the largest step for which a variety of simulations were stable and accurate. However, even when the small step of was implemented, the program ran faster by an order of magnitude (see Table 3).