A relaxed derivation of the SOR coefficient

We begin our discussion by considering a simpler iteration system,

$$x^{k+1}=Gx^k$$

Let C be the set of all sequences $\{x^k\}$ created by this iteration such that the sequence converges to $\{x^*\}$, i.e. we say

$$\lim_{k\rightarrow \infty}\{x^k\} = \{x^*\}$$

Define

$$\alpha=\lim_{k\rightarrow \infty}sup\vert x^k-x^*\vert^{1/k} \rightarrow 0 \leq \beta \leq 1$$ (3.11)

We shall take liberties and call this $\alpha$ the asymptotic convergence factor of the general iteration (see the references for a more sophisticated definition). For a given iteration, the sequence will converge at least as rapidly as the sequence $\alpha^k$ goes to zero.

From our previous definitions,

$$x[k+1]-x^*=H(x[k]-x^*) = H(H(x[k-1]-x^*)$$

$$\vert x^k - x^*\vert \leq \vert H^k\vert\vert x^0-x^*\vert \leq \... ...-x^*\vert \ni \lim_{k\rightarrow \infty}sup\vert x^k-x^*\vert^{1/k}\leq \rho(H)$$ (3.12)

Here we have used that $\vert H\vert\leq \rho(H)$. We have loosly demonstrated that

$$\alpha \leq \rho(H)$$

The spectral radius establishes the upper bound on the rate of convergence. We thus seek to maximize the spectral radius as a to minimize the convergence time. To begin to do so, we wish to roughly demonstrate a rather elegant lemma that says that eigenvalues of the Jacobi iteration matrix occur in pairs (except if there is a zero eigenvalue), which we shall denote $\pm \mu_i$

Suppose we take some random value $\lambda$ and want to find $det(I\lambda - C)$, where C has eigenvalues $a_i$. C can be diagonalized with its eigenvalues via $C=TDT^{-1}$, whereby

$$det(\lambda I -C)=\det(T\lambda I T^{-1} -T(Diag)T^{-1})=det(T)det(\lambda I - Diag)det(T^{-1})$$

Thus, supposing the existance of p zero valued eigenvalues (p may be zero itself):

$$det(\lambda I - C) = det(\lambda I - Diag) = \lambda^p \prod^r_{i=1}(\lambda - a_i)$$

For a more thorough discussion of what is to follow, the reader is referred to the references, especially on 2-Cyclic matrices (such as the one we are concerned with). Now $det(\lambda I - C)$ may be written as

$$det(\lambda I - C)=\left\vert\begin{array}{cc}\lambda I_1 & -C_1 \\ -C_2 & \lambda I_2\end{array}\right\vert$$

And we note that

$$(-1)^n = \left\vert\left(\begin{array}{cc}-I_1&0\\ 0&I_2\end{... ...t( \begin{array}{cc}I_1&0\\ 0&-I_2\end{array}\right)\right\vert$$

We can thus write that

$$det(\lambda I -C) = (-1)^n\left\vert\begin{array}{cc}-I_1&0\\ 0&I_2\end{array}\right\... ...ray}\right\vert \left\vert\begin{array}{cc}I_1&0\\ 0&-I_2\end{array}\right\vert$$

$$= (-1)^n\left\vert\begin{array}{cc}-\lambda I_1 & -C_1 \\ -C_2 & -\lambda I_2\end{array}\right\vert$$

$$= det(\lambda I +C)$$ (3.13)

$$det(\lambda I - C)=det(\lambda I + C) = \lambda ^p \prod^{r}_{i=1}(\lambda + a_i) \Rightarrow \forall a_i > 0 \ \exists \ a_s < 0 \ \ni a_s = -a_i$$ (3.14)

Thus (those with a strongly developed mathematical sense will be unhappy with my use of thus here) we have Romanovsky's Lemma:

$$det(\lambda I -C)=\lambda^p \prod_{i=1}^{r/2}(\lambda^2-\mu_i^2)$$

We can use this lemma to demonstrate another important detail about form of matrices we are considering (formally called regular splittings).

We wish to show that the eigenvalues of $\alpha D^{-1}L + \beta D^{-1}U$ are equal to those of $\alpha^{1/2}\beta^{1/2}D^{-1}(L+U)$. Let

$$\delta = \frac{\alpha^{1/2}}{\beta^{1/2}}$$

We need to demonstrate that

$$det(\lambda I - (\alpha D^{1}L + \beta D^{-1}U)) = det(\lambda I - \alpha^{1/2}\beta^{1/2}(\delta D^{-1}L + \delta^{-1}D^{-1}U))$$

$$= det(\lambda I - \alpha^{1/2}\beta^{1/2}(D^{-1}L + D^{-1}U)$$ (3.15)

That is to say, that the matrix $\eta = L+U$ is consistently ordered. Define

$$\eta(\delta)=\left(\begin{array}{cc}0 & \delta^{-1}D_1^{-1}C_1\\ \delta D_2^{-1}C_2 & 0 \end{array}\right)$$

And

$$S_{\delta}=\left(\begin{array}{cc}I_1 & 0 \\ 0 & \delta I_2 \end{... ...lta}=\left(\begin{array}{cc}I_1 & 0 \\ 0 & \delta^{-1} I_2 \end{array}\right)\$$

It is clear that

$$\eta(\delta)=S_{\delta}\eta(1)S_{\delta}^{-1}$$

Whereby in the spirit of previous arguments we have thus shown Equation 3.15 to be true. We are now ready to attack the problem head on. Recall that we defined

$$H = \left(D-wL\right)^{-1}\left((1-w)D+wU\right)$$ (3.16)

from this discussion we can reformulate this as follows:

$$det(\lambda I - H) = det((D-wL)^{-1})det(D-wL)det(\lambda I - \left(D-wL\right)^{-1}\left((1-w)D+wU\right))$$

$$= det((D-wL)^{-1})det((D-wL)\lambda I - (1-w)D+wU)$$

$$= det(D^{-1})det((D-wL)\lambda I - (1-w)D+wU)$$

$$= det((I-wD^{-1}L)\lambda I - (1-w)I - wD^{-1}U)$$

$$= det((\lambda + w -1)I -\lambda w D^{-1}L - wD^{-1}U)$$

$$= det((\lambda + w -1)I-\lambda^{1/2}wD^{-1}C) \ \ \ \ \alpha=w\lambda, \ \beta = w \$$

$$= \left(\lambda + w -1\right)^p \prod^{r/2}_{i=1}\left(\left(\lambda + w -1\right)^2I - \lambda w^2\mu_i^2\right) \$$ (3.17)

The eigenvalues of H may be found from

$$(\lambda + w -1)^2 = \lambda w^2 \mu_i^2$$

From this we get

$$(\sqrt\lambda)^2 \mp \sqrt\lambda w \mu_i + (w-1)=0 \ \rightarrow \sqrt\lambda = \frac{w\mu_i \pm \sqrt{w^2\mu_i^2-4(w-1)}}{2}$$

Since $\mu$ comes in $\pm$ pairs, we can minimize our convergence time by choosing w such that

$$w^2\mu_i^2 - 4(w-1)=0$$

We are bound by the maximal $\mu$, the spectral radius, so that

$$w^2\rho(D^{-1}C)^2-4w+4=0$$

Defining $D^{-1}C = J$ and solving for w, we find

$$w=\frac{2-2\sqrt{1-\rho(J)^2}}{\rho(J)^2}=\frac{2(1-1+\rho(J)^2)}{\rho(J)^2(1+\sqrt{1-\rho(J)}}=\frac{2}{1+\sqrt{1+\rho(J)^2}}$$ (3.18)

To find the spectral radius of J, we need find the eigenvectors/values of J. The Jacobi equation is given, again, by

$$4x_{ij}-\left(x_{i+1,j}+x_{i-1,j}+x_{i,j+1}+x_{i,j-1}\right)=b_{i,j}$$

J is formed by stripping off the diagonal and multiplying by $I/4$ Thus the eigenvalue equation gives

$$J\zeta = \nu \zeta \ \rightarrow \frac{1}{4}\left(\alpha_{i+1,j}+\alpha_{i-1,j}+\alpha_{i,j+1}+\alpha_{i,j-1}\right)=\nu \alpha_{i,j}$$ (3.19)

Introduce the ansatz that

$$\alpha_{i,j}^{k,l}=\sin(\frac{k\pi i}{N})\sin(\frac{l\pi j}{N})$$

where N is the dimensions of the grid (matrix, i.e. the number of rows). Plugging this into the equation above yields the equation

$$\nu_{k,l}=\frac{1}{2}\left(\cos(\frac{\pi k}{N})+\cos(\frac{\pi l}{N})\right)$$

Obviously the maximum possible value will be $\cos(\pi/N)$, this thus defining the spectral radius of concern, where by,

$$w=\frac{2}{1+\sqrt{1-\cos^2(\pi/N)}}=\frac{2}{1+\sin(\pi/N)}\approx\frac{2}{1+\frac{\pi}{N}}$$ (3.20)