o_O

April 12, 2008

Iterative Solvers in a Nutshell

Filed under: Uncategorized — secr @ pm

Basic concepts of iterative algorithms for elliptic PDEs can be compressed to the following remarks:

1. Continuous and discrete equations

(a) Given a PDE $Lu = v$ with Dirichlet boundary conditions, $L$ an linear elliptic differential operator, we wish to find the solution $u$ .

(b) The solution can be written formally as $u = L^{-1} v$ . $L^{-1}$ is a linear elliptic integral operator, given by $L^{-1} v(x) = \int G(y-x) \, v(y) \, dy$ , the convolution of the right hand side with some Green’s function of $L$ .

(c) $L$ and $v$ can be discretized (somehow) as $A$ (a matrix) and $b$ (a vector). The goal is to find the discretized version $x$ of $u$ , given by $A x = b$ .

(d) The formal solution of the discretized problem is $x = A^{-1} b$ . But inverting $A$ directly would be to expensive, that’s where the iterative solver springs into action.

2. The iterative solver

(a) The goal is to build an iteration $x_{i+1} = f(x_i, b)$ that converges against the solution $x = A^{-1} b$ .

(b) While $i$ increases, we wish that the error in $x = x_i$ ,

$e_i = x_i - x$ ,

decreases and finally reaches 0. We’ll never know the value of $e_i$ because $x$ is unknown. But that doesn’t matter as long as the error in $Ax_i = b$ ,

$r_i = b - A x_i$

is known. The name of $e_i$ is simply “the error” while $r_i$ is called “the residual”.

(c) The relation between the error and the residual is $- A e_i = r_i$ , i.e. $e_i$ scales with $r_i$ . The “anisotropic” scaling-factor is $A^{-1}$ .

(d) The relation between the exact solution and the residual is $x = x_i + A^{-1} r_i$ . If we replace $A^{-1}$ with an approximation $B$ , we’ll get an approximation of $x$ :

$x_{i+1} = x_i + B r_i$ .

I would frame this equation now if I could, because it defines $f(x_i, b)$ , the iterator that’ll hopefully converge against $x$ .

(e) To be sure that it does, $||e_i||$ has to approach 0 for $i \rightarrow \infty$ . The related proof is always done by induction, i.e. $e_{i+1}$ has to be related to $e_i$ . Note that our iteration $x_{i+1} = x_i + B r_i$ can be rewritten as $x_{i+1} = x_i - B A e_i$ (can be simply done with (c)), so that $e_{i+1} = x_{i+1} - x = x_i - B A e_i - x$ . So eventually

$e_{i+1} = (I - B A) e_i$ .

To conclude whether $x_i$ approaches $x$ , we have to make sure that $||e_i||$ approaches 0. The key is the following relation between $||e_{i+1}||$ and the spectral radius of $I - BA$ :

$||e_{i+1}|| \leq \rho(I-BA) ||e_i||$ .

You’ll see it if you expand $e_{i+1} = (I - B A) e_i$ in the Eigenbase of $I-BA$ . The conclusion is:

If $\rho(I - BA) < 1$ then convergence.

April 10, 2008

Heat Equation & Finite Differences

Filed under: Uncategorized — secr @ pm

The 2D-steady-state-heat-equation (hail to the hyphen) is

$\partial_x^2 u + \partial_y^2 u = f.$

Introducing a rectangular grid (step-width $h$ and $N$ points in each direction) on a finite domain (which need not to be rectangular) and substituting the continuous functions and operators with discretized versions

$\partial_x^2 u \rightarrow$ $\partial_x^{-} \partial_x^{+} u = (u_{m+1,l} - 2 u_{m,l} + u_{m-1,l})/h^2$

gives (provided the special case $N = 4$ and zero boundary conditions)

- Convergence ( $v_{ml} \rightarrow u(m,l)$ for $h \rightarrow 0$ ) can be shown (even though it’s a bit tricky, the proof makes use of the discrete version of the maximum-principle).

- The matrix belongs to the family of block-matrices, i.e. it’s built by (N-1)^2 blocks of quadratic sub-matrices of size N-1 (the “bandwith”), so the total number of elements is $(N-1)^4$ . Each block is either diagonal or tridiagonal. This tridiagonal-block-matrix shape is characteristic for all elliptic problems.

- Solving the system with the wrong algorithm can be highly expensive: Gauss-elimination costs $O(n^3)$ and bandsolvers $O(n^2)$ (with $n^2 = N$ ), while the price of the most sophisticated algorithms is $O(n \log n)$ (Buneman’s Total Reduction) or even $O(n)$ (adaptive multilevel Schwartz + Krylov).The problem with Total Reduction is that it can’t be applied to all elliptic problems. Good general-purpose solvers are iterative solvers, i.e. Jacobi, Gauss-Seidel, gradient method (all $O(n^2)$ ), SOR, conjugated gradient ( $O(n^{2/3})$ ) method plus some more sophisticated methods (adaptive multilevel Schwartz + Krylov, $O(n)$ ).

- In d dimensions: The bandwith of the matrix-blocks is given by $w=(N-1)^{(d-1)}$ , or approximately $n/N$ . The multilevel-method wins in all dimensions, but note that in 1D, bandsolvers aren’t worse.

full    Gauss elim.   bandsolver  Jacobi, GS, gradient    SOR, CG     fast solvers

d = 1        3            1                3                 2             1

d = 2        6            4                4                 3             2

d = 3        9            7                5                 4             3

d is the number of dimensions, the numbers # describe the asymptotic (i.e. large N) costs to reach a given accuracy as $O(N^\sharp)$ . In general, bandsolver’s cost is $O(n w^2) = O(N^d w^2)$ . With w = n/N it’s $O(N^{3d-2})$ , which is cheap for d=1, but expensive for d > 1.

o_O

April 12, 2008

Iterative Solvers in a Nutshell

April 10, 2008

Heat Equation & Finite Differences

Recent Posts

Pages

Follow “o_O”