Jacobi- and SOR method: visualization

Directly to the input form

• Jacobi's method and the SOR-(successive overrelaxation) method are two of the best known iterative solvers for systems of linear equations. Both methods transform a linear system of equations Ax=b by means of a so called ''splitting'' in a fixpoint equation x=Gx+g, where G and g may depend on a parameter ω which must be chosen properly. The fixed point equation is subsequently solved by direct iteration
  x^k+1 = Gx^k+g.
  It is easy to show that this iteration converges from arbitrary x⁰ if and only if ρ(G), the so called ''spectral radius'', i.e. the maximum of the absolute values of the eigenvalues of G, is less than one. This is obtained under quite restrictive conditions on A only. For Jacobi the iteration reads
  x^k+1 = x^k + D^-1(b-A*x^k)
  and for SOR
  (D-ω*L) x^k+1 = ((1-ω)*D + ω*U) x^k + ω*b
  where A is splitted into its strict lower triange -L, its diagonal D and its strict upper triangle -U. Computationally more convenient for SOR is the equivalent formulation
  x^k+1_i = x^k_i + (ω/A_i,i)(b- ∑_{j=1 to i-1} A_i,jx^k+1_j -∑_{j=i to n} A_i,jx^k_j), i=1,2,...,n
  The upper index denotes the iteration and the lower index the component of the vector. SOR with ω=1 is better known under the name ''Gauss-Seidel''. Jacobi and SOR require that no diagonal element of A vanishes. By induction over the column number it is easy to see that every invertible matrix can be permuted by rows such that at least this property is satisfied. But sufficient convergence criteria are restrictive: best known is the fact that SOR converges for symmetric positive definite matrices for 0 < ω < 2 whereas Jacobi doesn't necessarily converge in this case. Both methods converge for ''irreducibly diagonally dominant L-matrices'' for the parameter choice 0 < ω <= &omega₀, &omega₀ > 1 unknown, but unfortunately in this case the spectral radius decreases monotonically so that ω=1 is the best safe choice. A matrix is called an ''irreducibly diagonally dominant L-matrix'' if the following properties hold:
  
  • A_i,i > 0 , A_i,j <= 0 for i not equal j
    
  • A_i,i >= ∑_{j=1 to i-1,i+1 to n} |A_i,j| for all i
    with at least once a strict inequality
    
  • A cannot be transformed to upper block triangular form by a similarity transformation with a permutation matrix.
  Such a matrix is always invertible and has a positive inverse, it is a so called ''M-matrix''.
• Only under even more restrictive conditions, namely
  
  • The Jacobi method converges and its iteration matrix has only real eigenvalues
  • In the splitting A=-L+D-U the matrix
    D^-1(α*L+(1/α)*U) has eigenvalues independent on α (also formulated as ''the matrix is consistently ordered'')
  one knows that an unique optimal overrelaxation parameter ω exists which speeds up convergence considerable. It has the value
  ω_best= 2/ (1+ (1 - μ²)^1/2 )
  with μ the spectral radius of the Jacobi method (and in this case also μ² the spectral radius for Gauss-Seidel). In practice one doesn't know this μ and hence must also first approximate this.
• These methods are applied only if direct solvers cannot be used for reasons of memory or time limitations and if the solution is required to low precision only.
• The rate of convergence of Jacobi's method is in general worse than that of SOR with properly chosen parameter. In principle it is the spectral radius of the iteration matrix which determines this convergence rate, but the norm in which this applies is matrix dependent and often looks quite strange. Hence don't wonder if sometimes it has the appearance as if SOR would diverge although finally it comes down to the solution.
• For both methods the convergence rate depends on the condition number cond(A)=||A||||A^-1|| of the matrix A and you cannot hope for more than
  1 - const/cond(A)
  Only in the special case discussed above (consistently ordered matrix) SOR gives a rate
  1 - const/sqrt(cond(A)).
  
• In order to visualize the methods we use here the dimension n=2. In this case consistent ordering and irreducibility is satisfied if both nondiagonal elements are nonzero.

Input

• We use here n=2, hence the 6 data of the system are entered directly.
• You have the choice between the two methods and can specify the relaxation parameter for SOR. This must be real and in the open interval ]0,2[.
• The coefficients of the system are entered rowwise, the matrix elements followed by the right hand side. You must follow the conventions of FORTRAN. In case of a very bad conditioning of your matrix, more precisely if
  max{|A_i,j|}/|det(A)| > 100000
  this visualization makes no sense and you will get a corresponding message. But even if the conditioning is worse than 100 then you will see only one line and little if anything of the iteration sequence due to the limitations of the graphics resolution, but there is the printed table.
• you also specify the initial guess. For more visible effects this should be not a ''good'' one.
• You specify the termination criterion: if the euclidean length of the residual vector is less than ε with 1.E-8 <= ε <= 1.E-3 then we terminate.
• You also specify the maximum number of iterative steps allowed.
• Since the methods don't converge always, we need a criterion for detecting divergence. For this we use a ''blow up'' in the solution. This blow up is defined by leaving of a predefined box:
  ||x-x^*||_∞ <= (1+factor*||x⁰-x^*||_∞)
  Here x^* is the exact solution. You specify factor. In case of ω > ω_opt, then you must choose factor > 1. For ω <= 1 factor = 1 is sufficient. Important: 1.0 <= factor <= 1.E3 !

Output

• The linear system Ax=b is printed.
• You get a table containing the two components of x and the euclidean length of the residual r=b-Ax for all steps performed.
• A plot which visualizes the iteration: the two equations of the system are shown as lines in R² and of course their intersection is the solution sought. The points of the iteration sequence are shown, connected with arrows. Be aware that for Jacobi both corrections are applied simultaneously, hence you will see the arrows not path following the axis as for SOR.

Questions ?!

• When you get convergence, when divergence?
• n=2 is a quite special case: you might force convergence of Gauss-Seidel by a simple trick. Which one?
• How works the relaxation parameter ?
• In which case you must choose a large value for factor? What does that mean for the vectornorm in which convergence will be monotonic? (If there is convergence, then there is always a vectornorm, in which convergence is monotonic).
• Can you verify the dependence of the convergence rate from the condition number of the matrix?

To the input form