Numerical Linear AlgebraNumerical Linear AlgebraProbably the simplest kind of problem.Occurs in many contexts, often as part of larger problem.Symbolic manipulation packages can do linear algebra "analytically" (e.g. Mathematica, Maple).Numerical methods needed when:Number of equations very largeCoefficients all numericalLinear SystemsLinear Systems     =      = ⋮ ⋮     = Write linear system as:This system has n unknowns and m equations.If n = m, system is closed.If any equation is a linear combination of any others, equations are degenerate and system is singular.**see Singular Value Decomposition (SVD), NRiC 2.6.Numerical ConstraintsNumerical ConstraintsNumerical methods also have problems when:1) Equations are degenerate "within round-off error".2) Accumulated round-off errors swamp solution (magnitude of a's and x's varies wildly).For n,m < 50, single precision usually OK.For n,m < 200, double precision usually OK.For 200 < n,m < few thousand, solutions possible only for sparse systems (lots of a's zero).Matrix FormMatrix FormWrite system in matrix form:where: ==⋯  ⋯  ⋮ ⋮ ⋮ ⋯ ColumnsRowsMatrix Data RepresentationMatrix Data RepresentationRecall, C stores data in row-major form:a11, a12, ..., a1n; a21, a22, ..., a2n; ...; am1, am2, ..., amnIf using "pointer to array of pointers to rows" scheme in C, can reference entire rows by first index, e.g. 3rd row = a[2].Recall in C array indices start at zero!!FORTRAN stores data in column-major form:a11, a21, ..., am1; a12, a22, ..., am2; ...; a1n, a2n, ..., amnNote on Numerical Recipes in CNote on Numerical Recipes in CThe canned routines in NRiC make use of special functions defined in nrutil.c (header nrutil.h).In particular, arrays and matrices are allocated dynamically with indices starting at 1, not 0.If you want to interface with the NRiC routines, but prefer the C array index convention, pass arrays by subtracting 1 from the pointer address (i.e. pass p-1 instead of p) and pass matrices by using the functions convert_matrix() and free_convert_matrix() in nrutil.c (see NRiC 1.2 for more information).Tasks of Linear AlgebraTasks of Linear AlgebraWe will consider the following tasks:1) Solve Ax = b, given A and b.2) Solve Axi = bi for multiple bi's.3) Calculate A-1, where A-1A = I, the identity matrix.4) Calculate determinant of A, det(A).Large packages of routines available for these tasks, e.g. LINPACK, LAPACK (public domain); IMSL, NAG libraries (commercial).We will look at methods assuming n = m.The Augmented MatrixThe Augmented MatrixThe equation Ax = b can be generalized to a form better suited to efficient manipulation:The system can be solved by performing operations on the augmented matrix.The xi's are placeholders that can be omitted until the end of the computation. ∣ =⋯  ∣ ⋯  ∣ ⋮ ⋮ ⋮ ∣ ⋮  ⋯ ∣ Elementary Row OperationsElementary Row OperationsThe following row operations can be performed on an augmented matrix without changing the solution of the underlying system of equations:I. Interchange two rows.II. Multiply a row by a nonzero real number.III. Add a multiple of one row to another row.The idea is to apply these operations in sequence until the system of equations is trivially solved.The Generalized Matrix EquationThe Generalized Matrix EquationConsider the generalized linear matrix equation:Its solution simultaneously solves the linear sets:Ax1 = b1, Ax2 = b2, Ax3 = b3, and AY = I, where the xi's and bi's are column vectors.∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ =∣ ∣ ∣    ∣ ∣ ∣    ∣ ∣ ∣    ∣ ∣ ∣    Gauss-Jordan EliminationGauss-Jordan EliminationGJE uses one or more elementary row operations to reduce matrix A to the identity matrix.The RHS of the generalized equation becomes the solution set and Y becomes A-1.Disadvantages:1) Requires all bi's to be stored and manipulated at same time ⇒ memory hog.2) Don't always need A-1.Other methods more efficient, but good backup.Gauss-Jordan Elimination: ProcedureGauss-Jordan Elimination: ProcedureStart with simple augmented matrix as example:Divide first row (a1|b1) by first element a11.Subtract ai1 (a1|b1) from all other rows:Continue process for 2nd row, etc.∣ ∣ ∣  / / ∣ /  −/  − /  ∣ −/  −/  − /  ∣ −/ Row a1|b1Pivot rowFirst column of identity matrixGJE Procedure, Cont'dGJE Procedure, Cont'dProblem occurs if leading diagonal element ever becomes zero.Also, procedure is numerically unstable!Solution: use "pivoting" - rearrange remaining rows (partial pivoting) or rows & columns (full pivoting - requires permutation!) so largest coefficient is in diagonal position.Best to "normalize" equations (implicit pivoting).Gaussian Elimination with Gaussian Elimination with BacksubstitutionBacksubstitutionIf, during GJE,