### Linear Systems of Equations Iterative Methods

```Systems of Linear Equations
Iterative Methods
Ax  b
Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften
ETH Hönggerberg / HCI F128 – Zürich
E-Mail: [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
1
Iterative Methods
 The main idea behind iterative methods is that we can
reformulate the problem in the following way
Ax  S  (A  S) x  b
S x  b  (S  A ) x
 Now we can proceed iteratively
( k  1)
(k )
Sx
 b  (S  A ) x
 This is advantageous if S has a structure that makes it
particularly easy to factorize or solve (e.g. diagonal,
tridiagonal, triangular) and if S can be considered constant
 Even if b or A (and S) change slightly, the solution will be
similar and thus the iteration will converge quickly if the
previous solution is used as a starting point
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
2
Convergence of Iterative Methods
 Let us rewrite the convergence loop as
x
( k  1)
 Mx
(k )
c

 M  S 1  S  A 

1
c

S
b

 As we have seen with the convergence loop of a predictor
corrector ODE-solver scheme, such a loop will only
converge if the spectral radius of the iteration matrix is
smaller than 1, i.e.
1
  M     I  S A    m ax  1
 This is the case if S and A are similar in some sense:
S
1
 A
1

1
S A I
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
3
Solution Procedure
 If we look at the iteration equation
Sx
( k  1)
 b  (S  A ) x
(k )
we can formulate it in a more convenient way by treating
the right hand side as a constant vector
( k  1)
Sx
c
 Since we can choose S to have a structure that is easy to
solve, there is no need to calculate S-1 or do Gauss
elimination to solve this equation system; We can use more
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
4
Jacobi’s Method
 The Jacobi method chooses S to be a diagonal matrix
Dx
( k  1)
 b  (D  A ) x
(k )
 diag ( A )
D
 The procedure is therefore given by
D
 a1,1

a
 2 ,1
A

 a N 1,1
a
 N ,1
a1,2
a 2 ,2
a1,3
a 2 ,3
a N 1, N  2
a N ,N 2
a N 1, N 1
a N , N 1
a1, N 

a2,N



a N 1, N 
a N , N 
(0)
xi

( k  1)
xi
bi
a i ,i

bi
a i ,i

1
a i ,i

(k )
ai, j x j
ji
 This method is guaranteed to converge if A is diagonally
dominant, i.e.
a 
a
i ,i

i, j
ji
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
5
The Gauss-Seidel Method
 In this case, S is chosen to be lower triangular
Lx
( k  1)
 b  (L  A ) x
(k )
 tril ( A )
L
 The procedure is therefore given by
L
 a1,1

a
 2 ,1
A

 a N 1,1
a
 N ,1
a1,2
a1,3
a 2 ,2
a 2 ,3
a N 1, N  2
a N 1, N 1
a N ,N 2
a N , N 1
a1, N 

a2,N



a N 1, N 
a N , N 
x
( k  1)
i

1 
(k )
( k  1)

 bi   a i , j x j   a i , j x j 
a i ,i 
ji
ji

 This method is guaranteed to converge if A is symmetric
positive-definite or diagonally dominant
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
6
Another way of Implementation
 At every iteration, a linear system of equations is solved:
Sx
( k  1)
 b  (S  A ) x
(k )
 This is equivalent to
1
S  A 
M
S
c
S b
x
1
( k  1)
 Mx
(k )
c
 It is therefore possible to compute M and c and once in the
beginning and then do simple matrix multiplications and
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
7
 We say that two non-zero vectors are conjugate (wrt A) if
u Av  0
 If A is symmetric and positive definite, this corresponds to
an inner product of u and v
T
 If we now find a series of N mutually conjugate vectors p,
the solution of the linear system Ax = b must be contained
in the space that these vectors span, i.e.
N
x

i
pi
i 1
 With the coefficients
T
N
b

i 1
i
A pi

k 
pk b
T
pk A pk
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
8
 If we find the vectors p one after the other, we can proceed
iteratively
 The main idea for the transformation into an iterative
method is that the unique minimizer of the cost function
f (x) 
1
x Ax  x b
T
T
2
is the same as the solution of the equation Ax = b
 The negative gradient of f in point xk, which denotes the
rk  b  A x k
 However to ensure that our next step be conjugate to the
previous one, we will have to correct this search direction
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
9
Stability and Convergence of the CG-Method
 The CG-Method is only stable if A is symmetric and
positive definite; However, if we have an asymmetric
matrix, we can use the equivalent normal equations
A Ax  A b
T
T
Cx
 z
where C = ATA is symmetric positive definite for any nonsingular matrix A
 The convergence speed of the CG-Method is determined
by the condition number of A; The larger it is, the slower
the method converges and unfortunately
 A A   
T
2
A 
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
10
 We can improve the convergence speed by
preconditioning the equation system, i.e. replace the
system by an equivalent system
Ax  b

M
1
1
Ax  M b
so that κ(M-1A) < κ(A)
 Some choices of M are
 Jacobi-preconditioning where M = D; D is the diagonal matrix of A
 Incomplete Cholesky factorization where M = KKT and KKT ≈ A
 SSOR-preconditioning where for ω = (0,...,2)
 1
 1 
M 
D

L

 D 
2 
  
1
1
 1

D

L




T
where L is the strictly lower triangular matrix of A
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
11
Algorithm for the PCG-Method
1. Choose a starting point x0 and calculate
r0  b  A x 0 ;
z 0  M \ r0 ;
2. Proceed from xk to xk+1 using pk as the direction
T
k 
rk z k
T
pk A pk
x k 1  x k   k p k
rk  1  rk   k A p k
3. Correct the search direction
z k  1  M \ rk  1
T
k 
rk  1 z k  1
T
rk z k
Iterate 2 and 3 until
norm(rk+1) or norm(Axk – b)
is sufficiently small, or
5*length(A) steps have
been taken
p k 1  z k 1   k p k
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
12
Exercise: 2-D Heat Transport
 Consider a spherical catalyst particle where a chemical
reaction takes place
 Assumptions:




The reaction rate is independent of concentration and temperature
Thermal diffusivity is independent on temperature
No convective heat transport
Perfect heat sink at the boundary, i.e. T = 0
λ
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
13
Exercise (Continued)
 The heat transfer can be described as a PDE:
  2Q  2Q 
2





Qk
2
2 
t
y 
 x
Q
Q  x  y  r ,t  0
2
2
2
Q  x, y, t  0   0
where  2 is the Laplace operator, k is the heat produced by
the reaction and r is the particle radius
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
14
Exercise (Continued)
 Substituting the temperature:
Q   c PT
T
t
  T 
2
k
 cP
 At steady state we get
k
  cP
  T
2
 This equation can be solved by using a discretized Laplace
operator:
DT 
k
  cP
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
15
Assignment
1. Implement the PCG-algorithm as a funciton (see slide 12)
2. Use it to solve the discrete steady state Laplace equation
for the catalyst particle, discretizing with N = 50 points
 Assume
k
  cP
1
 Create a disk shaped grid using: G = numgrid('D',N);
 Check the grid with spy(G)
 Use D = delsq(G); to create a negative 5th-order discrete
Laplace operator, take a look at the operator again with spy(D)
 Create the right hand side using b = ones(size(D,1),1);
 Choose an initial guess, e.g. x0 = zeros(size(D,1),1);
3. Solve the system without preconditioning
 Set M = eye(size(D));
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
16
Assignment (continued)
4. Plot the solution using
 U = G;
U(G>0) = full(x(G(G>0)));
clabel(contour(U))
5. Try out different preconditioning methods:
 Jacobi preconditioning, i.e. M = diag(diag(D));
 SSOR preconditioning with ω = 1.5; Use L = tril(D,-1); to
get the strictly lower triagonal part and E = diag(diag(D)); to
get the diagonal matrix of D
 Incomplete cholesky factorization; This can be done using
M = C*C';
with
C = ichol(D); OR
C = cholinc(D,'0');
 How many iterations are needed in all cases?
Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods
17
```