The Jacobi Method

We first divide the matrix $\bss{A}$ into 3 parts

$$\bss{A} = \bss{D} + \bss{L} + \bss{U}$$ (3.13)

where $\bss{D}$ is a diagonal matrix (i.e. $D_{ij} = 0$ for $i\not= j$) and $\bss{L}$ and $\bss{U}$ are strict lower and upper triangular matrices respectively (i.e. $L_{ii} = U_{ii} = 0$ for all $i$). We now write the Jacobi or Gauss-Jacobi iterative procedure to solve our system of equations as

$$\bss{X}^{n+1} = \bss{D}^{-1}\left[\bss{B} - \left(\bss{L} + \bss{U}\right)\bss{X}^n\right]$$ (3.14)

where the superscripts ${}^n$ refer to the iteration number. Note that in practice this procedure requires the storage of the diagonal matrix, $\bss{D}$, and a function to multiply the vector, $\bss{X}^n$ by $\bss{L} + \bss{U}$. This algorithm resembles the iterative solution of hyperbolic or parabolic partial differential equations, and can be analysed in the same spirit. In particular care must be taken that the method is stable. Simple C code to implement this for a 1D Poisson's equation is given below.

      int i;
      const int N = ??;  // incomplete code
      double xa[N], xb[N], b[N];
      while ( ... ) // incomplete code
      {
         for ( i = 0; i < N; i++ )
            xa[i] = (b[i] - xb[i-1] - xb[i+1]) * 0.5;
         for ( i = 0; i < N; i++ )
            xb[i] = xa[i];
      }

Note that 2 arrays are required for X and that the matrices, $\bss{D}$, $\bss{L}$ and $\bss{U}$ don't appear explicitely.