Summary of Relaxation Methods

Eqn. 13 is the starting point for defining four common relaxation strategies. The two factors which distinguish these methods are the numerical value of $\omega$ and whether the updated values of the function are used once they are obtained. The direct use of Eqn. 13 is the weighted Jacobi method. If $\omega$ is set to one, then it is called the Jacobi method. In this case, the updated function is just a local average of the neighboring function values plus a term from the fixed total charge on the grid point i.

The second relaxation type uses the updated value of the function once obtained, i.e. $\phi_{i-1}^t$ is replaced by $\phi_{i-1}^{t+1}$.This may seem like a minor change, but it actually accelerates the convergence considerably. If a similar spectral analysis is made to that above, one finds that $\omega$ values between one and two can be used, which implies a larger `time step' size. If $\omega = 1$, it is called Gauss-Seidel iteration, and for $\omega$ values different than one, it is Successive Over-Relaxation, or SOR. Generally SOR is the most efficient means of relaxing on a given grid, with $\omega = 1.5$ a suitable choice for one dimensional iterations. The best values of $\omega$ may differ in higher dimensions depending on how it is defined. A few trial and error attempts will show the optimal relaxation weighting parameter; if it is too large, the relaxation becomes unstable. Typically the optimal value is slightly less than the maximum one allowed for stability. SOR methods are not necessarily strictly variational, since the updated function is used for the next grid point iteration. However, in practice the action generally goes downhill as with the strict steepest descents and line minimization algorithms.

Before relaxation steps are implemented, one must first make a choice of the boundary conditions (fixed or periodic), and the initial guess for the $\phi$ function. The boundaries should be set according to the physics of the problem. If a finite system is being examined, the most common case will be to set the function to zero on the boundaries, as long as there is some screening to make the potential short ranged. For periodic or `wrap-around' boundary conditions (PBC), the value on one side is set exactly the same as that on the opposite side. We note that there is no unique solution for the PBC case; any constant can be added to the $\phi$ field with no physical consequence (as long as there is charge balance, which is required for stability).

Once the boundary conditions are set, an initial guess is made [19], and a relaxation strategy is chosen, the iterations begin. One simply `sweeps' through the lattice in some specific order, thus updating all N grid points. After a given sweep, how can we monitor whether our new approximation is good enough? If we write the matrix equation to be solved as Lu=f where L is the Laplacian, u the current potential approximation, and f minus the grid charge, then we have a natural estimate of how close we are to the grid solution. Define r=f-Lu, which is termed the residual; then, when |r| = 0, the problem has been solved exactly on the lattice (there still remain errors between the lattice approximation and the `true' continuum differential equation). The residual should decrease with each iteration in order to have a convergent algorithm, and the iterations can be terminated when the magnitude of the residual has reached a certain limiting cutoff based on a physical criterion, say convergence of the system total energy in relation to the thermal energy k_B T (k_B is Boltzmann's constant and T is the absolute temperature).

The relaxation equation (Eqn. (13)) and thus the four relaxation strategies can be easily extended to higher dimensions. Take the two dimensional example. In this case, the analog of Eq. (13) becomes:

$$\phi_{i,j}^{t+1} = (1-\omega) \phi_{i,j}^{t} +\frac{\omega}{... ...i_{i+1,j}^t+\phi_{i,j+1}^t + \frac{4\pi}{\epsilon} \rho_{i,j}).$$ (17)

Here $\rho_{i,j}$ is the charge density per unit transverse distance (i.e., the 3-d charge density times h²). Then the lattice sweeps are over all grid points (i,j). The four types of relaxation can be defined precisely as for the one dimensional case depending on the specific update scheme and value of the relaxation parameter $\omega$ [20].