1Numerical Methods in Linear Numerical Methods in Linear Algebra, Part TwoAlgebra, Part TwoLarry CarettoMechanical Engineering 501ASeminar in Engineering Seminar in Engineering AnalysisAnalysisSeptember 16, 20042Outline• Review last class• Pivoting strategies to reduce round-off error in solutions• Gauss-Jordan for matrix inverses• The LU method• Numerical analysis software3Review last class• Finite representation of numbers in binary computer gives round-off error• Machine epsilon is measure of precision of floating point representation–1 + ε = 1 below “machine epsilon”– epsilon is usually 1.19x10-7 for single and 2.22x10-16for double precision• Round-off error growth in calculations called error propagation4Review Error Propagation• Addition and subtraction errorsyxyxyxrel~~−=εxyyxxyrel~~−=εyxyxrel~~εεε+=• Multiplication and division errors• Same result for both()()yxyxyxyxyxyxyxrel~~~~±+≈±+=±±−±≡εεεεε5Review Matrix Vector Norm• Both Ax and x are matrices• Can use any matrix norm to compute ||A|| xAxAxmax=• Choosing infinity norm as vector norm gives row sum• Choosing one norm as vector norm gives column sum∑==niijja11maxA∑=∞=njijia1maxA6Review Problem Condition• In ill-conditioned problems small relative changes in data cause large relative changes in results• Matrix condition number κ(A) = ||A|| ||A-1|| of 100 or above indicate ill-conditioningAδAAxδx)(κ≤bδbAxδx)(κ≤brAxxx)(~κ≤−27Review Gauss Elimination• Use each row from row 1 to row n-1 as the “pivot” row– Work on each row below the pivot row• Multiply pivot row by arow,pivot/apivot,pivot• Subtract result from row r to make arow,pivot= 0• Operation requires subtraction for each column of A right of pivot column and for b– Repeat for each row below pivot• Repeat for rows 1 to n-1 as pivot rows• Use back substitution for x values8Review Round-off Error Example• Same set of equations solved with original order and reverse order⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡0002.11300003.01110002.111300003.0xyyx⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡−=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡− 9999.019997.20110003.00002.110002.199990300003.0xyyx• Original order solution less accurate by about four significant figures9Review Round-off Conclusions• Showed that order of equations was important factor in round-off error• Problems caused by small pivot elements (diagonal element on pivot row)• Found large loss of significant figures with original order but no error when order was reversed• Want to use this idea in algorithms for reducing round-off error10Zero and Small Pivots• Gauss elimination may give zero for the pivot element in one row, but swapping equations can give a nonzero pivot• Always check for a small number not zero• Use scaled equations so maximum absolute element in each row is one• Find row with maximum (scaled) element in the pivot column and then swap11Pivoting Code for Homeworkp = k ‘k is row to be pivot rowbig = Abs(a(k, k) / s(k))For i = k + 1 To NIf (Abs(a(i, k) / s(i)) > big) Thenbig = Abs(a(i, k) / s(i))p = iEnd IfNext ipivot = big < tol ‘singular matrix flag‘Swap rows p and k to get maximum on pivot12Modified Code for Homework‘Make sure there is not a zero pivotFor i = k + 1 To Np = ibig = Abs(a(i, k) / s(k))If (big > tol) ThenExit ForEnd IfNext ipivot = big <= tol ‘singular matrix flag‘Swap rows p and k to get nonzero on pivot313Finding Inverse Matrices•For B = A-1, AB = I, gives⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1000010000100001321333323122322211131211321333323122322211131211LLMOMMMMOMMMLLLLLLLLMOMMMMOMMMLLLLLLLLMOMMMMOMMMLLLLLLnnnnnnnnnnnnnnnnbbbbbbbbbbbbbbbbaaaaaaaaaaaaaaaa• This requires n solutions of n equations in n unknowns (one solution for each b column) 14Finding Inverse Matrices II• E. g., for the second column of B we solve⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡00102322212321333323122322211131211MMMMLLMOMMMMOMMMLLLLLLnnnnnnnnnbbbbaaaaaaaaaaaaaaaa• We can solve for all n columns of the inverse matrix at one time15Finding Inverse Matrices III• Gauss-Jordan subtracts pivot row from all rows above and below to give⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡232221223222121000010000100001nnbbbbββββMMMMLLMOMMMMOMMMLLLLLL• Diagonal form gives solutions by inspection: b12= β12, b22= β22, b32= β32, …, bn2= βn2,16Gauss-Jordan• Start with augmented A matrix⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1000010000100001321333323122322211131211LLMOMMMMOMMMLLLLLLLLMOMMMMOMMMLLLLLLnnnnnnnnaaaaaaaaaaaaaaaa• Gauss elimination with diagonals set to 1 and pivot row subtracted from all rows 17Gauss-Jordan Result• Augmented matrix at end of process⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡nnnnnnnnbbbbbbbbbbbbbbbbLLMOMMMMOMMMLLLLLLLLMOMMMMOMMMLLLLLL3213333231223222111312111000010000100001• Inverse is read from augmented (right-hand) part of matrix 18Gauss-Jordan PseudocodeLoop over all rows as pivots (1 to N)Find the row with the maximum (scaled) element in the pivot column and swap it with the pivot current pivot rowDivide all elements the pivot row by a[pivot,pivot]Loop over all rows, subtracting the pivot row times a[row,pivot] from each row (1 to N except pivot)419Gauss-Jordan Example• Solve the set of equations on the right)(14837)(13923)(342642321321321iiixxxiixxxixxx=++=++−−=−−• Divide equation (i) by a11= 2)(14837)(13923)(171321321321321iiixxxiixxxixxx=++=++−−=−−• Subtract –3 times (i) from equation (ii) and 7 times (i) from (iii) to replace (ii) and (iii)20Gauss-Jordan Example II()[]()[]()[]()[])17(313)13(39)2(32133321−−−=−−−+−−−+−−−xxx()[]()[]()[]()[])17(714)13(78)2(73177321−−=−−+−−+−xxx• Result from first set of operations)(13399170)(383040)(171321321321321iiixxxiixxxixxx=++−=−−−=−−• Divide equation