Computational Mechanics

Iterative Solver and an overview of Domain Decomposition Methods

Iterative Solver

First, let's talk about iterative solver. Some basic iterative techniques and their properties will be introuced.

Outline

Stationary iterative methods (relaxation, smoothing)

• Richardson
• Jacobi
• Gauss-siedel
• SOR

Projection methods

• Conjugate Gradient
• GMRES and variants

Basic iterative methods❗️

The first iterative methods we will discuss are the basic iterative
methods. Basic iterative methods only use information of the previous
iteration. They are also called stationary method, or relaxation
methods or smoothing methods

👆Most basic iterative methods can be constructed as a splitting of A :

$A=M−N$

M is usually chosen as a matrix that is easy to invert, whereas N is the rest.

The original linear system $Ax = B$ plug in the decomposed expression, we can get $(M - N) x = b$. Reconstruct this formula, we obtain $M x=N x+b$. We introduce iteration metrics $k$ and $k + 1$ to represent the current step and next step solution estimates. So we get the iterative formula:

$M x^{k+1}=N x^k+b$

The above iterative process can be rewritten in an equivalent form to make it easier to understand and implement:

$x^{k+1}=M^{-1}\left(N x^k+b\right)$

where $N = M - A$, then we obtain:

$x^{k+1}=x^k+M^{-1}\left(b-A x^k\right)$

$x^{k+1}=x^k+M^{-1} r^k$

where $\mathbf{r}^k=\mathbf{b}-\mathbf{A} \mathbf{x}^k$ is the residual for iterate $k$. And the Matrix $\mathbf{M}$ is a preconditioner.

🌟Richardson method

Let's take $\mathbf{M}=\mathbf{I}$ as our preconditionner. Then $\mathbf{N}=\mathbf{I}-\mathbf{A}$.

From basic iterative method, we know $\mathbf{x}^{k+1}=\mathbf{x}^k+\mathbf{M}^{-1} \mathbf{r}^k$. Then we use $I$ to replace $M$:

$\begin{aligned} \mathbf{x}^{k+1} & =\mathbf{x}^k+\mathbf{r}^k \\ & =\mathbf{x}^k+\mathbf{b}-\mathbf{A} \mathbf{x}^k \\ & =\mathbf{b}+(\mathbf{I}-\mathbf{A}) \mathbf{x}^k \end{aligned}$

Initial guess $\mathbf{x}^0=\mathbf{0}$, then we have:

$\begin{aligned} & \mathbf{x}^1=\mathbf{b} \\ & \mathbf{x}^2=\mathbf{b}+(\mathbf{I}-\mathbf{A}) \mathbf{x}^1=\mathbf{b}+(\mathbf{I}-\mathbf{A}) \mathbf{b} \\ & \mathbf{x}^3=\mathbf{b}+(\mathbf{I}-\mathbf{A}) \mathbf{x}^2=\mathbf{b}+(\mathbf{I}-\mathbf{A}) \mathbf{b}+(\mathbf{I}-\mathbf{A})^2 \mathbf{b} \end{aligned}$

We rewrite these items as:

$\mathbf{x}^{k+1}=\sum_{i=0}^k(\mathbf{I}-\mathbf{A})^i \mathbf{b}$

This is actually the expression of our linear system $A x = b$, so we can deduce $\sum_{i=0}^k(\mathbf{I}-\mathbf{A})^i$ should represent $A^{-1}$.

🔥We define inverse matrix as $A A^{-1} = I$, so let's have a look at $\sum_{i=0}^{\infty}(\mathbf{I}-\mathbf{A})^i$, if it exist. Let $\mathbf{B}=\mathbf{I}-\mathbf{A}$

we construct $A A^{-1}$ by $\mathbf{A} \sum_{i=0}^k(\mathbf{I}-\mathbf{A})^i$, then we introduce $\mathbf{B} = \mathbf{I} - \mathbf{A}$ to this formula:

$\sum_{i=0}^k\mathbf{B}^i (\mathbf{I} - \mathbf{B})$

To evaluate this we need to use an emprical method:

For $k = 2$:

$\begin{aligned} & \sum_{i=0}^2(I-A)^i A \\ = & \left(I+B+B^2\right)(I-B) \\ = & I+B+B^2-B-B^2-B^3 \\ = & I-B^3 \end{aligned}$

Then we deduce:

$\sum_{i=0}^k\mathbf{B}^i (\mathbf{I} - \mathbf{B}) = \mathbf{I} - \mathbf{B}^{n+1}$

It is equivalent to:

$\mathbf{A} \sum_{i=0}^k(\mathbf{I}-\mathbf{A})^i = \mathbf{I} - (\mathbf{I} - \mathbf{A})^{n+1}$

So if Richardson converges, when $n$ goes to infinity, $(\mathbf{I} - \mathbf{A})^{n+1}$ should go to zero. That leads to our discussion about the convergence condition of Richardson.

🌈The method will converge if $(\mathbf{I}-\mathbf{A})^i$ converges. Let $\lambda$ be any eigen value of $\mathbf{A},(\mathbf{I}-\mathbf{A})^i$ will converge if $|1-\lambda|<1$. In other word, we need $\rho(\mathbf{I}-\mathbf{A})=\rho\left(\mathbf{M}^{-1} \mathbf{N}\right)<1$, the spectral radius of the iterative matrix must be less than 1. If the eigen values are real, this means $0<\lambda<2$

The iterative process continues until a stopping criterion is reached, such as the norm of residuals is less than a given threshold, or the number of iterations reaches a given maximum.

The Richardson method is a simple iterative method, and its convergence mainly depends on the condition number of the matrix A. Generally, if the condition number of A is smaller, the Richardson method will converge faster. However, if the condition number of A is large, then convergence may be slow or may not even converge. In practical applications, convergence may be improved by choosing a different preprocessor or using other iterative methods

🌟Jacobi Method