Classical iterative algorithms

by Yi Zhang

For this topic, many excellent textbooks are available. Most of the modern numerical PDE solvers rely on one or another iterative algorithm taking care of the linear system given by discretization ( I know very little of direct methods). In this scenario, one has linear system

Ax=b

and tries to come up with some iterative scheme

x^{(k+1)}=Tx^{(k)}+c

with

\lim_{k\to\infty}{x^{(k)}}=x^{\ast}

meaning the convergence to the “real” solution, or, to be exact, the solution by direct solver if we disregard the round-off error. By this setup, one can see the solution is the fixed point of operator T. So the convergence of the algorithm is equivalent to two proposals: the spectral radius \rho(A)<1, and there exists operator norm induced by some norm of the linear space to which x belongs that this norm of T is less than 1. The convergence speed is measured by average rate of convergence of A (for m iterations) defined by

R(A^m)=-\log[\|A^m\|^{1/m}]=-\log[\|A^m\|]/m

which comes from the fact that \lim_{k\to\infty}(\|T^k\|^{1/k})=\rho(T). In order to put the original problem into this fixed point setup, A is splitted:

A=M-N\longrightarrow x^{(k+1)}=M^{-1}Nx^k+M^{-1}b

If A is splitted according to entries’ positions as A=L+D+U, one has the classical algorithms:

Jacobi: M=D,N=-(L+U)

Gauss-Seidel: M=D+L,N=-U

SOR: x_i^{(k+1)}=\omega \bar{x}_i^{(k+1)}+(1-\omega)x_i^{(k)}

The last one is in the component wise form in which \bar{x}^{(k+1)} is the Gauss-Seidel result, and the relaxation parameter is between zero and two (even though the name “overrelaxation” comes when it is greater than one). Each of those algorithms has direct(but different in implementation) analogical block version. All these methods are stationary, i.e. in each iteration T is unchanged. If this is not the case, we have unstationary methods, in which after each iteration information is re-gathered in order to determine a new T for next iteration. When A is a SPD matrix, a large family of unstationary methods can be formed by the solving a minimization problem(minimal residual method (MINRES)). Classical CG method I wrote about before belongs to this family. A modern member of this family is GMRES by Saad and Schultz in 1986.

Advertisements