Lecture notes on Numerical Linear Algebra

numerical linear algebra in signals systems and control and numerical linear algebra for high-performance computers pdf free download
LaylaChadwick Profile Pic
LaylaChadwick,New Zealand,Teacher
Published Date:15-07-2017
Your Website URL(Optional)
Comment
NUMERICAL LINEAR ALGEBRA Rolf Rannacher Institute of Applied Mathematics Heidelberg University Lecture Notes WS 2013/2014 February 7, 2014ii Address of Author: Institute of Applied Mathematics Heidelberg University Im Neuenheimer Feld 293/294 D-69120 Heidelberg, Germany rannacheriwr.uni-heidelberg.de http://www.numerik.uni-hd.deContents 0 Introduction 1 0.1 Basic notation of Linear Algebra and Analysis. . . . . . . . . . . . . . . . . . . . 1 0.2 Linear algebraic systems and eigenvalue problems . . . . . . . . . . . . . . . . . . 2 0.3 Numerical approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0.4 Applications and origin of problems . . . . . . . . . . . . . . . . . . . . . . . . . 4 0.4.1 Gaussian equalization calculus . . . . . . . . . . . . . . . . . . . . . . . . 4 0.4.2 Discretization of elliptic PDEs . . . . . . . . . . . . . . . . . . . . . . . . 5 0.4.3 Hydrodynamic stability analysis . . . . . . . . . . . . . . . . . . . . . . . 9 1 Matrix Analysis: Linear Systems and Eigenvalue Problems 13 n 1.1 The normed Euclidean space K . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.1 Vector norms and scalar products . . . . . . . . . . . . . . . . . . . . . . 13 1.1.2 Linear mappings and matrices . . . . . . . . . . . . . . . . . . . . . . . . 21 1.1.3 Non-quadratic linear systems . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.1.4 Eigenvalues and eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.1.5 Similarity transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.1.6 Matrix analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.2 Spectra and pseudo-spectra of matrices . . . . . . . . . . . . . . . . . . . . . . . 35 1.2.1 Stability of dynamical systems . . . . . . . . . . . . . . . . . . . . . . . . 35 1.2.2 Pseudospectrum of a matrix . . . . . . . . . . . . . . . . . . . . . . . . . . 38 1.3 Perturbation theory and conditioning. . . . . . . . . . . . . . . . . . . . . . . . . 43 1.3.1 Conditioning of linear algebraic systems . . . . . . . . . . . . . . . . . . . 43 1.3.2 Conditioning of eigenvalue problems . . . . . . . . . . . . . . . . . . . . . 45 1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2 Direct Solution Methods 53 2.1 Gaussian elimination, LR and Cholesky decomposition . . . . . . . . . . . . . . . 53 2.1.1 Gaussian elimination and LR decomposition . . . . . . . . . . . . . . . . . 53 2.1.2 Accuracy improvement by final iteration . . . . . . . . . . . . . . . . . . . 61 2.1.3 Inverse computation and the Gauß-Jordan algorithm . . . . . . . . . . . . 63 2.2 Special matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 2.2.1 Band matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 iiiiv CONTENTS 2.2.2 Diagonally dominant matrices . . . . . . . . . . . . . . . . . . . . . . . . . 69 2.2.3 Positive definite matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 2.3 Irregular linear systems and QR decomposition . . . . . . . . . . . . . . . . . . . 72 2.3.1 Householder algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 2.4 Singular value decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 2.5 “Direct” determination of eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . 83 2.5.1 Reduction methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 2.5.2 Hyman’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 2.5.3 Sturm’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 2.6 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3 Iterative Methods for Linear Algebraic Systems 93 3.1 Fixed-point iteration and defect correction . . . . . . . . . . . . . . . . . . . . . . 93 3.1.1 Stopping criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.1.2 Construction of iterative methods . . . . . . . . . . . . . . . . . . . . . . 98 3.1.3 Jacobi- and Gauß-Seidel methods . . . . . . . . . . . . . . . . . . . . . . . 101 3.2 Acceleration methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.2.1 SOR method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.2.2 Chebyshev acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 3.3 Descent methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 3.3.1 Gradient method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 3.3.2 Conjugate gradient method (CG method) . . . . . . . . . . . . . . . . . . 121 3.3.3 Generalized CG methods and Krylov space methods . . . . . . . . . . . . 126 3.3.4 Preconditioning (PCG methods) . . . . . . . . . . . . . . . . . . . . . . . 128 3.4 A model problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 3.5 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 4 Iterative Methods for Eigenvalue Problems 141 4.1 Methods for the partial eigenvalue problem . . . . . . . . . . . . . . . . . . . . . 141 4.1.1 The “Power Method” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 4.1.2 The “Inverse Iteration” . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 4.2 Methods for the full eigenvalue problem . . . . . . . . . . . . . . . . . . . . . . . 146 4.2.1 The LR and QR method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 4.2.2 Computation of the singular value decomposition . . . . . . . . . . . . . . 151CONTENTS v 4.3 Krylov space methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 4.3.1 Lanczos and Arnoldi method . . . . . . . . . . . . . . . . . . . . . . . . . 154 4.3.2 Computation of the pseudospectrum . . . . . . . . . . . . . . . . . . . . . 159 4.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5 Multigrid Methods 173 5.1 Multigrid methods for linear systems . . . . . . . . . . . . . . . . . . . . . . . . . 173 5.1.1 Multigrid methods in the finite element context . . . . . . . . . . . . . . . 174 5.1.2 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 5.2 Multigrid methods for eigenvalue problems (a short review) . . . . . . . . . . . . 187 5.2.1 Direct multigrid approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 5.2.2 Accelerated Arnoldi and Lanczos method . . . . . . . . . . . . . . . . . . 189 5.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 Bibliography 194 Index 198vi CONTENTS0 Introduction Subject of this course are numerical algorithms for solving problems in Linear Algebra, such as linear algebraic systems and corresponding matrix eigenvalue problems. The emphasis is on iterative methods suitable for large-scale problems arising, e. g., in the discretization of partial differential equations and in network problems. 0.1 Basic notation of Linear Algebra and Analysis Atfirst,weintroducesomestandardnotation inthecontext of(finitedimensional)vector spaces of functions and their derivatives. Let K denote the field of real or complex numbers R or C, n respectively. Accordingly, for n ∈ N, let K denote the n-dimensional vector space of n- tuples x = (x ,...,x ) with components x ∈K, i = 1,...,n. For these addition and scalar 1 n i multiplication are defined by: x+y :=(x +y ,...,x +y ), αx :=(αx ,...,αx ), α∈K. 1 1 n n 1 n n The elements x ∈ K are, depending on the suitable interpretation, addressed as “points” or “vectors” . Here, one may imagine x as the end point of a vector attached at the origin of the 1 chosen cartesian coordinate systemandthe components x asits “coordinates” with respectto i this coordinate system. In general, we considervectors as “column vectors”. Within the “vector T calculus” its row version is written as (x ,...,x ) . The zero vector (0,...,0) may also briefly 1 n be written as 0. Usually, we prefer this coordinate-oriented notation over a coordinate-free 1 k n notation because of its greater clearness. A set of vectors a ,...,a in K is called “linearly independent” if k X i ca =0, c ∈K ⇒ c =0, i=1,...,k. i i i i=1 n Such a set of k =n linearly independent vectors is called a “basis” of K , which spans all of n n K , i. e., each element x∈K can (uniquely) be written as a linear combination of the form n X i x = ca, c ∈K. i i i=1 n Each (finite dimensional) vector space, such as K , possesses a basis. The special “carte- 1 n i sian basis” e ,...,e is formed by the “cartesian unit vectors” e := (δ ,...,δ ), δ = 1 1i ni ii and δ = 0, for i = 6 j, being the usual Kronecker symbol. The elements of this basis ij are mutually orthonormal, i. e., with respect to the euclidian scalar product, there holds P n j i j i n×n (e,e ) := e e = δ . “Matrices” A ∈ K are two-dimensional square arrays of 2 ij k=1 k k n numbers from K written in the form A =(a ) , where the first index, i, refers to the row ij i,j=1 and the second one, j, to the column (counted fromthe left uppercorner of the array) at which the element a is positioned. Usually, matrices are square arrays, but in some situations also ij rectangular matrices mayoccur. Thesetof(square) matrices formsavector space withaddition 1 Ren´e Descartes (1596-1650): French mathematician and philosopher (“cogito ergo sum”); worked in the Netherlands and later in Stockholm; first to recognize the close relation between geometry and arithmetic and founded analytic geometry. 12 Introduction and scalar multiplication defined in the natural elementwise sense, n n n A=(a ) , B =(b ) , c∈K ⇒ cA+B =(ca +b ) . ij ij ij ij i,j=1 i,j=1 i,j=1 For matrices and vectors natural multiplications are defined by d d     X n X n n n×n Ax= a x ∈K , AB = a b ∈K . ik k ik kj i=1 i,j=1 k=1 k=1 d Matrices are used to represent linear mappings in K with respect to a given basis, mostly a T T n ¯ cartesian basis, ϕ(x) = Ax. By A = (a ) , we denote the conjugate “transpose” of a ij i,j=1 n n×n T n×n matrix A = (a ) ∈K with the elements a = a¯ . For matrices A,B ∈K there ij ji i,j=1 ij T T T T ¯ holds (AB) =B A . Matrices for which A=A are called “symmetric” in the case K=R and “hermitian” in the case K=C. n Further,wehavetodealwithscalarorvector-valued functions u =u(x)∈K forarguments n x∈K . For derivatives of differentiable functions, we use the notation 2 2 ∂u ∂ u ∂u ∂ u 2 2 ∂ u := , ∂ u:= , ... , ∂u:= , ∂ u := , ... , x i x ij 2 ∂x ∂ x ∂x ∂x∂x i i j and analogously also for higher-order derivatives. With the nabla operator ∇ the “gradient” of a scalar function and the “divergence” of a vector function are written as gradu = ∇u := d (∂ u,...,∂ u) and divu = ∇·u := ∂ u +...+∂ u , respectively. For a vector β ∈ R the 1 d 1 1 d d derivative in direction β is written as ∂ u :=β·∇u. Combination of gradient and divergence β 2 yields the “Laplacian -Operator” 2 2 ∇·∇u=Δu=∂ u+...+∂ u. 1 d m The symbol ∇ u denotes the “tensor” of all partial derivatives of order m of u, i. e., in two j 2 i dimensions u=u(x ,x ), ∇ u=(∂ ∂ u) . 1 2 i+j=2 1 2 0.2 Linear algebraic systems and eigenvalue problems Let A be an m×n-matrix and b an m-vector,     a ··· a b 11 1n 1     . . . m,n m     . . . A=(a ) = , b=(b ) = . jk . . j . j,k=1 j=1     a ··· a b m1 mn m 2 PierreSimonMarquisdeLaplace(1749-1827): Frenchmathematician andastronomer; prof. inParis; founded among other fields probability calculus.0.3 Numerical approaches 3 We seek an n-vector x=(x ) such that k k=1,...,n a x + a x +··· + a x = b 11 1 12 2 1n n 1 . . (0.2.1) . a x +a x +··· +a x = b m1 1 m2 2 mn n m or written in short as Ax = b. This is called a “linear system” (of equations). It is called “underdetermined” for m n, “quadratic” for m = n, and “overdetermined” for m n. The linear system is solvable if and only if rank(A) =rankA,b (rank(A) = number of linearly independent columns of A) with the composed matrix   a ··· a b 11 1n 1   . . .   . . . A,b = . . . .   a ··· a b m1 mn m In the “quadratic” case the solvability of the system (0.2.1) is equivalent to any one of the n×n following properties of the coefficient matrix A∈K : - Ax =0 implies x=0. - rank(A) =n. - det(A)6=0. - All eigenvalues of A are nonzero. n×n A number λ∈ C is called “eigenvalue” of the (quadratic) matrix A∈ K if there exists a n corresponding vector w∈K \0, called ”eigenvector”, such that Aw =λw. (0.2.2) Eigenvalues are just the zeros of the characteristic polynomial χ (z) := det(A−zI) of A, A so that by the fundamental theorem of Algebra each n×n-matrix has exactly n eigenvalues counted accordingly to their (algebraic) multiplicities. The corresponding eigenvectors span lin- n ear subspaces of K called “eigenspaces”. Eigenvalue problems play an important role in many problems from science and engineering, e. g., they represent energy levels in physical models (e. g., Schro¨dinger equation in Quantum Mechanics) or determine the stability or instability of dynamical systems (e. g., Navier-Stokes equations in hydrodynamics). 0.3 Numerical approaches We will mainly consider numerical methods for solving quadratic linear systems and associated eigenvalue problems. Theemphasis will beon medium- and large-scale problems, i. e., problems 4 9 of dimension n≈10 −10 , which at the upper end pose particularly hard requirements on the algorithms with respect to storage and work efficiency. Problems of that size usually involve matrices with special structure such as “band structure” and/or extreme “sparsity”, i. e., only very few matrix elements in each row are non-zero. Most of the classical methods, which have originally been designed for “full” but smaller matrices, cannot realistically be applied to such large problems. Therefore, modernmethodsextensively exploit theparticular sparsity structure ofthematrices. Thesemethodssplitintotwoclasses, “directmethods”and“iterative methods”.4 Introduction Definition 0.1: A “direct” method for the solution of a linear system Ax =b is an algorithm, which (neglecting round-off errors) delivers the exact solution x in finitely many arithmetic operations. “Gaussian elimination” isatypical example ofsucha “directmethod”. In contrast to t that an “iterative method” constructs a sequenceof approximate solutions x , which only in t∈N t the limit t→∞ converge to the exact solution, i. e., lim x =x. “Richardson iteration” or t→∞ more general fixed-point methods of similar kind are typical example of such “iterative methods”. In analyzing a direct method, we are mainly interested in the work count, i. e., the asymptotic number of arithmetic operations needed for achieving the final result depending on the problem 3 size, e. g., O(n ), while in an iterative method, we look at the work count needed for one iteration step and the number of iteration steps for reducing the initial error by a certain fixed −1 factor, e. g., 10 , or the asymptotic speed of convergence (“linear”, “quadratic, etc.). However, there is no sharp separation between the two classes of “direct” or “iterative” methodsasmanytheoretically“direct”methodsareactuallyusedin“iterative” forminpractice. Atypical methodofthistypeistheclassical “conjugate gradient (CG)method”ofHestenes and Stiefel, which in principle is a direct method (after n iteration steps) but is usually terminated like an iterative methods already after m≪ n steps. 0.4 Applications and origin of problems We present some applications, from which large linear algebra problems originate. This illus- trates how the various possible structures of matrices may look like. 0.4.1 Gaussian equalization calculus A classical application in Astronomy is the Gaussian equalization calculus (method of least 1 n 2 error-squares): For given functions u ,...,u and points (x ,y )∈R , j = 1,...,m, mn, j j a linear combination n X k u(x) = c u (x) k k=1 is to be determined such that the “mean deviation” m   X 1/2 2 Δ := u(x )−y 2 j j j=1 3 becomes minimal. (The so-called “Chebyshev equalization problem” in which the “maximal deviation” Δ := max u(x )−y isminimized posesmuch moresevere difficulties and ∞ j=1,...,m j j is therefore used only for smaller n.) For the solution of the Gaussian equalization problem, we set y :=(y ,...,y ), c≡(c ,...,c ), 1 m 1 n k k k 1 n a :=(u (x ),...,u (x )), k =1,...,n, A≡a ,...,a . 1 m 3 Pafnuty Lvovich Chebyshev (1821-1894): Russian mathematician; prof. in St. Petersburg; contributions to number theory, probability theory and especially to approximation theory; developed the general theory of orthogonal polynomials.0.4 Applications and origin of problems 5 Using this notation, now the quadratic functional m   X 1/2 2 F(c) = (Ac−y) j j=1 n is to be minimized with respect to c ∈ R . This is equivalent to solving the overdetermined linear system Ac = y in the sense of finding a vector c with minimal mean error-squares, i. e., with minimal “defect”. In case that rank(A) = n this “minimal-defect solution” c is determined by the so-called “normal equation” T T A Ac =A y, (0.4.3) a quadratic linear n×n-system with a positive definite (and hence regular) coefficient matrix T k k−1 A A. In the particular case of polynomial fitting, i. e., u (x) =x , the “optimal” solution n X k−1 u(x) = c x k k=1 is called “Gaussian equalization parabola” for the points (x ,y ), j =1,...,m. Because of the j j 4 regularity of the “Vandermondian determinant”   n−1 1 x ··· x 1 1   n−1 n   Y 1 x ··· x 2 2   det  = (x −x ) 6= 0, k j . . .  . . .  . . . j,k=1,jk   n−1 1 x ··· x n n for mutually distinct points x there holds rank(A) = n, i. e., the equalization parabola is j uniquely determined. 0.4.2 Discretization of elliptic PDEs The numerical solution of partial differential equations requires first an appropriate “discretiza- tion” of the differential operator, e. g., by a “difference approximation” of the derivatives. Consider, for example, the so-called “first boundary value problem of the Laplacian operator”, Lu:=−Δu=f in Ω, u=g on ∂Ω, (0.4.4) n posed on a domain Ω ⊂ R with boundary ∂Ω. Here, for a given right-hand side function f =f(x ,x ) and boundary function g =g(x ,x ), assumed to be continuous, an on Ω twice 1 2 1 2 ¯ differentiable and on Ω continuous solution function u = u(x ,x ) is to be determined such 1 2 ¯ that (0.4.4) holds. The region Ω, e. g., the unit square, is covered by an equidistant cartesian mesh Ω with“meshboundary” ∂Ω andthemeshpoints P ∈Ω maybenumberedrow-wise. h h h 4 Alexandre-Thophile Vandermonde (1735-1796): French mathematician; gifted musician, came late to math- ematics and published here only four papers (nevertheless member of the Academy of Sciences in Paris); contri- butions to theory of determinants and combinatorial problem (curiously enough the determinant called after him does not appear explicitly in his papers).6 Introduction 1 2 3 4 1 5 6 7 8 h= mesh width m+1 9 10 11 12 2 n=m “inner” mesh points 13 14 16 15 h h Figure 1: Finite difference mesh At “interior” mesh points P ∈ Ω the differential operators in x - and x -direction are 1 2 h approximated by second-order central difference quotients, which act on mesh functions u (P). h This results in “difference equations” of the form X L u (P) := σ(P,Q)u (Q) =f (P), P ∈Ω , (0.4.5) h h h h h Q∈N(P) u (P) =g (P), P ∈∂Ω , (0.4.6) h h h with certain mesh neighborhoods N(P)⊂Ω∪∂Ω ofpoints P ∈Ω and approximations f (·) h h h to f and g (·) to g. Wesetthecoefficients σ(P,Q) :=0 forpoints Q6∈N(P). Theconsidered h difference operator based on second-order central difference quotients for approximating second derivatives is called “5-point difference operator” since it uses 5 points (Accordingly, its three- dimensional analogue is called “7-point difference operator”). Then, for P ∈Ω there holds h X X σ(P,Q)u (Q) =f (P)− σ(P,Q)g (Q). (0.4.7) h h h Q∈Ω Q∈∂Ω h h For any numbering of the mesh points in Ω and ∂Ω , Ω =P, i=1,...,n, ∂Ω =P, i = h h h i h i n+1,...,n+m, we obtain a quadratic linear system for the vector of approximate mesh values N U =(U ) , U :=u (P ). i i i h i=1 AU =F, (0.4.8) n n with A =(a ) , F =(b ) , where ij j i,j=1 j=1 n+m X a :=σ(P,P ), b :=f (P )− σ(P ,P )g (P ). ij i j j h j j i h i i=n+1 In the considered special case of the unit square and row-wise numbering of the interior mesh points Ω the 5-point difference approximation of the Laplacian yields the following sparse h 2 matrix of dimension n=m :0.4 Applications and origin of problems 7       B −I 4 −1 m m               −I B −I −1 4 −1 m m m 1     A = n B = m,     m . . 2 . .     h . . −I B −1 4 m m           . . . .   . . . . . . . . where I is them×m-unit matrix. The matrix A is a very sparse band matrix with half-band m width m, symmetric and (irreducibly) diagonally dominant. This implies that it is regular and 3 positivedefinite. Inthreedimensionsthecorrespondingmatrixhasdimension n=m andhalf- 2 band width m and shares all the other mentioned properties of its two-dimensional analogue. 4 7 In practice, n ≫ 10 up to n ≈ 10 in three dimensions. If problem (0.4.4) is only part of a larger mathematical model involving complex domains and several physical quantities such as (in chemically reacting flow models) velocity, pressure, density, temperature and chemical 7 9 species, the dimension of the complete system may quickly reach up to n≈10 −10 . To estimate a realistic size of the algebraic problem oriented by the needs of a practical application, we consider the above model problem (Poisson equation on the unit square) with an adjusted right-hand side and boundary function such that the exact solution is given by u(x,y) =sin(πx)sin(πy), 2 −Δu=2π u =:f in Ω, u=0 on ∂Ω. (0.4.9) For this setting the error analysis of the difference approximation yields the a priori estimate 1 2 2 2 maxu−u ≈ d M (u)h ≈8h , (0.4.10) h 4 Ω 24 Ω h 4 4 where M (u) = max ∇ u≈π (see the lecture notes Rannacher 3). In order to guarantee ¯ 4 Ω −3 −2 4 a relative error below TOL = 10 , we have to choose h≈ 10 corresponding to n≈ 10 in 6 two and n≈ 10 in three dimension. The concrete structure of the matrix A depends on the numbering of mesh points used: (i) Row-wise numbering: The lexicographical ordering of mesh points, (x,y ) ≤ (x ,y ) i j p q −1 for j≤q or j =q, i≤p, leads to a band matrix with band width 2m+1≈h . The sparsity within the band would be largely reduced by Gaussian elimination (so-called “fill-in”) 21 22 23 24 25 16 17 18 19 20 11 12 13 14 15 6 7 8 9 10 1 2 3 4 5 Figure 2: Lexicographical ordering of mesh points8 Introduction (ii) Diagonal numbering: The successive numbering diagonally to the cartesian coordinate directions leads to a band matrix with less band volume. This results in less fill-in within Gaussian elimination. 11 16 20 23 25 7 12 17 21 24 4 8 13 18 22 5 9 14 19 2 1 3 6 10 15 Figure 3: Diagonal mesh-point numbering (iii) Checker-board numbering: The staggered row-wise and column-wise numbering leads −1 to a 2×2-block matrix with diagonal main blocks and band width 2m+1≈h . 11 24 12 25 13 21 9 22 10 23 6 19 7 20 8 16 4 17 5 18 1 14 2 15 3 Figure 4: Checkerboard mesh-point numbering 5 For large linear systems of dimension n≫ 10 direct methods such as Gaussian elimination are difficult to realize since they are generally very storage and work demanding. For a matrix 6 2 of dimension n = 10 and band width m = 10 Gaussian elimination requires already about 8 10 storage places. This is particularly undesirable if also the band is sparse as in the above example with at most 5 non-zero elements per row. In this case those iterative methods are more attractive, in which essentially only matrix-vector multiplications occur with matrices of similar sparsity pattern as that of A. As illustrative examples, we consider simple fixed-point iterations for solving a linear system Ax=b with a regular n×n-coefficient matrix. The linear system is rewritten in the form n X a x + a x =b , j =1,...,n. jj j jk k j k=1 k6=j0.4 Applications and origin of problems 9 If a 6=0, this is equivalent to jj n n o X 1 x = b − a x , j =1,...,n. j j jk k a jj k=1 k6=j t n Then, the so-called “Jacobi method” generates iterates x ∈ R , t = 1,2,..., by successively solving n n o X 1 t t−1 x = b − a x , j =1,...,n. (0.4.11) j jk j k a jj k=1 k6=j t t Whencomputing x theprecedingcomponents x ,rj, arealreadyknown. Hence, inorderto j r accelerate the convergence of the method, one may use this new information in the computation t 5 of x . This idea leads to the “Gauß-Seidel method”: j n o X X 1 t t t−1 x = b − a x − a x , j =1,...,n. (0.4.12) j jk jk j k k a jj kj kj The Gauß-Seidel method has the same arithmetic complexity as the Jacobi method but under certain conditions (satisfied in the above model situation) it converges twice as fast. However, though very simple and maximal storage economical both methods, Jacobi as well as Gauß- Seidel, are by far too slow in practical applications. Much more efficient iterative methods are the Krylov-space methods. The best known examples are the classical “conjugate gradient method” (“CG method”) of Hestenes and Stiefel forsolving linear systems with positive definite matrices and the “Arnoldi method” for solving corresponding eigenvalue problems. Iterative methodswithminimalcomplexitycanbeconstructedusingmulti-scaleconcepts(e.g.,geometric or algebraic “multigrid methods”). The latter type of methods will be discussed below. 0.4.3 Hydrodynamic stability analysis Another origin of large-scale eigenvalue problems is hydrodynamicstability analysis. Let vˆ,pˆ be a solution (the “base flow”) of the stationary Navier-Stokes equation −νΔvˆ+vˆ·∇vˆ+∇pˆ=0, ∇·vˆ=0, in Ω, (0.4.13) in vˆ =0, vˆ =v , ν∂ vˆ−pˆn =P, vˆ =q, Γ Γ n Γ Γ out rigid in Q where vˆ is the velocity vector field of the flow, pˆ its hydrostatic pressure, ν the kinematic viscosity (for normalized density ρ ≡ 1), and q the control pressure. The flow is driven by a in prescribed flow velocity v at the Dirichlet (inflow) boundary (at the left end), a prescribed meanpressure P attheNeumann(outflow)boundary(attherightend)andthecontrolpressure at the control boundary Γ . The(artificial) “free outflow” (also called “do nothing”) boundary Q condition in (0.4.13) has proven successful especially in modeling pipe flow since it is satisfied by Poiseuille flow (see Heywood et al. 39). 5 Philipp Ludwig von Seidel (1821-1896): German mathematician; prof. in Munich; contributions to analysis (method of least error-squares) and celestial mechanics and astronomy.10 Introduction . Γ Q Γ S Γ in out Γ Q Figure 5: Configuration of the flow control problem. Figure 5 shows the configuration of a channel flow around an obstacle controlled by pressure prescription at Γ , and Figure 6 the computational mesh and streamline plots of two flows for Q different Reynolds numbers and control values, one stable and one unstable. Figure 6: Computational mesh(top), uncontrolled stable (middle)and controlled unstable(bot- tom) stationary channel flow around an obstacle. For deciding whether these base flows are stable or unstable, within the usual linearized stability analysis, oneinvestigates thefollowingeigenvalue problemcorrespondingtotheNavier- Stokes operator linearized about the considered base flow: −νΔv+vˆ·∇v+v·∇vˆ+∇q =λv, ∇·v =0, in Ω, (0.4.14) v =0, ν∂ v−qn =0. Γ ∪Γ n Γ out rigid in Fromthelocation oftheeigenvaluesinthecomplexplane,onecandrawthefollowingconclusion: Ifaneigenvalue λ∈C of(0.4.14) has Re λ0, thebaseflowisunstable,otherwiseitissaidto be “linearly stable”. This means that the solution of the linearized nonstationary perturbation problem ∂ w−νΔw+vˆ·∇w+w·∇vˆ+∇q =0, ∇·w =0, in Ω, t (0.4.15) w =0, ν∂ w−qn =0 Γ ∪Γ n Γ rigid in out0.4 Applications and origin of problems 11 corresponding to an initial perturbation w =w satisfies a bound 0 t=0 supkw(t)k≤Akw k, (0.4.16) 0 t≥0 withsomeconstant A≥1. Afterdiscretizationtheeigenvalueproblem(0.4.14)infunctionspace is translated into an essentially nonsymmetric algebraic eigenvalue problem, which is usually of 5 6 high dimension n≈10 −10 . Therefore its solution can only be achieved by iterative methods. However, “linear stability” doesnot guarantee full“nonlinear stability” dueto effects caused by the “non-normality” of the operator governing problem (0.4.14), which may cause the con- stant A to become large. This is related to the possible “deficiency” (discrepancy of geometric and algebraic multiplicity) or a large “pseudospectrum” (range of large resolvent norm) of the criticaleigenvalue. Thiseffectiscommonlyacceptedasexplanationofthediscrepancyinthesta- bility propertiesof simple base flows such as Couette flowand Poiseuille flow predicted by linear eigenvalue-based stability analysis and experimental observation (see, e.g., Trefethen&Embree 19 and the literature cited therein).12 Introduction1 Matrix Analysis: Linear Systems and Eigenvalue Problems In this chapter, we introduce the basic notation and facts about the normed real or complex n n×n vector spaces K of n-dimensional vectors and K of corresponding n×n-matrices. The n emphasis is on square matrices as special representations of linear mappings in K and their spectral properties. n 1.1 The normed Euclidean space K 1.1.1 Vector norms and scalar products We recall some basic topological properties of the finite dimensional “normed” (vector) space n K , where depending on the concrete situation K=R (real space) or K=C (complex space). n In the following each point x ∈ K is expressed by its canonical coordinate representation 1 n n x=(x , ...,x ) in terms of a (fixed) cartesian basis e ,...,e of K , 1 n n X i x= xe. i i=1 n Definition 1.1: Amapping k·k :K →R isa“(vector)norm”ifithasthefollowingproperties: n (N1) Definiteness: kxk≥0, kxk =0 ⇒ x=0, x∈K . n (N2) Homogeneity: kαxk =αkxk, α∈K, x∈K . n (N3) Triangle inequality: kx+yk≤kxk+kyk, x,y∈K . The notion of a “norm” can be defined on any vector space V over K, finite or infinite dimen- sional. The resulting pair V,k·k is called “normed space”. Remark 1.1: The propertykxk≥ 0 is a consequence of the other conditions. With (N2), we obtain 0 =k0k and then with (N3) and (N2) 0 =kx−xk≤kxk+k−xk = 2kxk. With the help of (N3) we obtain the useful inequality n kxk−kyk ≤kx−yk, x,y∈K . (1.1.1) Example 1.1: The standard example of a vector norm is the “euclidian norm” n   X 1/2 2 kxk := x . 2 i i=1 The first two norm properties, (N1) and (N2), are obvious, while the triangle inequality is a special case of the “Minkowski inequality” provided below in Lemma 1.4. Other examples of useful norms are the “maximum norm” (or “l norm”) and the “l norm” ∞ 1 n X kxk := max x, kxk := x. ∞ i 1 i i=1,...,n i=1 1314 Matrix Analysis: Linear Systems and Eigenvalue Problems The norm properties of k·k and k·k are immediate consequences of the corresponding ∞ 1 properties of the modulus function. Between l norm and maximum norm there are the so- 1 called “l norms” for 1p∞: p n   X 1/p p kxk := x . p i i=1 Again the first two norm properties, (N1) and (N2), are obvious and the triangle inequality is the Minkowski inequality provided in Lemma 1.4, below. ′ ′ n With the aid of a norm k·k the “distance” d(x,x) := kx−xk of two vectors in K is defined. This allows the definition of the usual topological terms “open”, “closed”, “compact”, n “diameter”, and “neighborhood” for point sets in K in analogy to the corresponding situation in K. We use the maximum norm k·k in the following discussion, but we will see later that ∞ n this is independent of the chosen norm. For any a∈K and r0, we use the ball n K (a) :=x∈K : kx−ak r r ∞ as standard neighborhoodof a with radius r. Thisneighborhoodis “open”since for each point c x∈K (a) there exists a neighborhood K (x)⊂K (a); accordingly the complement K (a) is r δ r r “closed”. The “closure” of K (a) is defined by K (a) :=K (a)∪∂K (a) with the “boundary” r r r r ∂K (a) of K (a). r r k n Definition 1.2: A sequence of vectors (x ) in K is called k∈N k - “bounded” if all its elements are contained in a ball K (0), i. e., kx k ≤R, k∈N, R ∞ k l - “Cauchy sequence” if for each ε∈R there is an N ∈N such that kx −xk ε, k,l≥N , + ε ∞ ε n k - “convergent” towards an x∈K if kx −xk →0 (k→∞). ∞ k k k For a convergent sequence (x ) , we also write lim x = x or x → x (k → ∞). k∈N k→∞ Geometrically thismeansthat anystandardneighborhood K (x) of x contains almost all (i. e., ε k all but finitely many) elements x . This notion of “convergence” is obviously equivalent to the componentwise convergence: k k kx −xk →0 (k→∞) ⇔ x →x (k→∞), i =1,...,n. ∞ i i n This allows the reduction of the convergence of sequences of vectors in K to that of sequences of numbers in K. As basic results, we obtainn-dimensional versions of the Cauchy criterion for convergence and the theorem of Bolzano-Weierstraß. Theorem 1.1 (Theorems of Cauchy and Bolzano-Weierstraß): n n (i) Each Cauchy sequence in K is convergent, i. e., the normed space (K ,k·k ) is complete ∞ (a so-called “Banach space” ). n (ii) Each bounded sequence in K contains a convergent subsequence. k n Proof. (i) For any Cauchy sequence (x ) , in view of x≤kxk , i=1,...,n, for x∈K , k∈N i ∞ k also the component sequences (x ) , i = 1,...,n, are Cauchy sequences in K and therefore k∈N i

Advise: Why You Wasting Money in Costly SEO Tools, Use World's Best Free SEO Tool Ubersuggest.