Iterative Approaches

An alternative method to solve a system of linear equations is iteratively through the use of a recursive relation, for example  

$${\bf B} \cdot {\bf x}^{n+1} ={\bf C} \cdot {\bf x}^n + {\bf b}$$ (244)

If the matrices ${\bf B}$ and ${\bf C}$ are selected such that ${\bf A} = {\bf B} - {\bf C}$, it is obvious from Eq. (8.6) that if it converges, the solution satisfies the original set of linear equations ${\bf A} \cdot {\bf x} = {\bf b}$. The conditions for convergence for this iterative algorithm are the same as those applicable, in general, for a fixed point algorithm for the solution of a nonlinear set of equations (see below), i.e., they are determined by the requirements of the contraction mapping theorem, which, in the present case, requires the norm of the Jacobian matrix of the transformation, $\vert\vert {\bf B^{-1} } \cdot {\bf C} \vert\vert < 1$. In addition, the matrices ${\bf B}$ and ${\bf C}$ are determined based on the computational load required to solve the recursive equation (8.6). Popular choices involve the Jacobi method (where ${\bf B}$ represents the diagonal component of ${\bf A}$ ) and the Gauss- Seidel method (where ${\bf B}$represents the lower triangular (including the diagonal) component of ${\bf A}$ ) both resulting in a recursive expression involving an order n² number of operations for a dense matrix ${\bf A}$. Thus for either one to succeed, less than n iterations should be required for convergence. It is easy to show that convergence in this case is guaranteed if the matrix ${\bf A}$ is diagonally dominant (i.e., the magnitude of any diagonal element is larger that the sum of absolute values of the rest of the elements belonging to the same row or column. Unfortunately, however, even in that case, a fast convergence is not always guaranteed. Indeed, if the matrix ${\bf A}$arises from the discretization of a differential equation, as is most often the case, it has been observed that as the size of the discretization increases, so does the number of iterations.

Several alternatives have been proposed. The Successive OverRelaxation method (SOR) uses a mixture of Jacobi and Gauss-Seidel methods with the relative weight chosen to be a fitted parameter. Another approach, and more elegant for complex problems, involves the use of a preconditioner: The original matrix equation is now premultiplied by a suitably chosen matrix so that to make the Jacobian matrix problem easier to solve. This is a practice that usually works with most of the solution approaches where an iterative solution of a nonlinear problem is sought. Of course, the best preconditioner would have involved the inverse of the matrix ${\bf A}$ itself, but this is exactly what one tries to solve, which means that one would have to settle with an approximation to the inverse. Several approximations have been proposed for the inverse, the best ones involving the first few steps of other iterative techniques. Finally, a technique that has become one of the most popular given its relatively straightforward implementation and theoretical justification is conjugate gradients. See [2] for more details.