Report

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] 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 direct (even analytical) approaches instead 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 ji This method is guaranteed to converge if A is diagonally dominant, i.e. a a i ,i i, j ji 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 ji ji 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 vector additions for the iteration Daniel Baur / Numerical Methods for Chemical Engineers / Iterative Methods 7 The Conjugate Gradient Method 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 The Conjugate Gradient Method (Continued) 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 steepest descent direction, reads 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 Preconditioned Conjugate Gradient 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 Qk 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