= 3 1.204 Lecture 20 Linear systems: Gaussian elimination Gaussian elimination LU decomposition Systems of Linear Equations 3x0+ x1 -2x2 2x2x0 +4x + 4x1 + 3x2+3x x0 -3x1 1 -2 2 4 3 2 4 3 1 -3 0 x0 x x11 x2 = = 5 =35 35 = -5 5 35 35 -5 A x = b 3 x 3 3 x 1 3 x 1 1Algorithm to Solve Linear System =aijCreate matrix x0 x1 x22 = b0 b1 b22 A x b 0 0 0 =a’ij x0 x1 x2 = Forward solve b’0 b’1 b’2 A’A’ xx b’b’ 0 0 0 x0Back solve x0 x1 x2 = b’0 b’1 b’2 x1 x2 A’ x b’ Gaussian Elimination: Forward Solve 3 1 -2 5 Q= 2 4 3 35Q 1 -3 0 -5 A b Form Q for convenience Do elementary row ops: Multiply rows Multiply rows Add/subtract rows Make column 0 have zeros below diagonal Pivot= 1/3 Pivot= 2/3 3 1 -2 0 10/3 13/3 00 -10/310/3 2/32/3 5 95/3 Row 1’= row 1 - (2/3) row 0 20/3 Row 2’=row 2 (1/3) row 0 Row 2=row 2 -(1/3) row 0-20/3 Make column 1 have zeros below diagonal 3 1 -2 5 Pivot= -1 0 10/3 13/3 95/3 Row 2’’= row 2’ + 1 * row 1 0 0 15/3 75/3 2o o00 Gaussian Elimination: Back Solve 3 1 -2 5 0 10/3 13/3 95/30 10/3 13/3 95/3 (15/3)x2=(75/3) x2 = 50 0 15/3 75/3 3 1 -2 5 0 10/3 13/3 95/30 0 15/3 75/3 (10/3)x1 + (13/3)*5= (95/3) x1 = 3 0 10/3 13/3 95/30 0 15/3 75/3 3x0 + 1*3 - 2*5 = 5 x0 = 43 1 -2 5 A Complication 0 1 -2 5 22 44 33 35 35 Row 1’= row 1 -(2/0) row 0(/ )o 1 -3 0 -5 Exchange rows: put largest pivot element in row: 2 4 3 35 00 11 -22 55 1 -3 0 -5 Do this as we process each column. If there is no nonzero element in a column, matrix is not full rank. 3Gaussian Elimination public class Gauss { public static double[] gaussian(double[][] a, double[] b) { int n = a length; int n = a.length; // Number of unknowns // Number of unknowns double[][] q = new double[n][n + 1]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) // Form q matrix q[i][j]= a[i][j]; q[i][n]= b[i]; } forward_solve(q); // Do Gaussian elimination back_solve(q); // Perform back substitution double[] x= new double[n]; // Extract column n of q, for (int i = 0; i < n; i++) // which contains the solution x x[i]= q[i][n]; return x; } Forward Solve private static void forward_solve(double[][] q) { int n = q.length; int m= q[0].length; for (int i = 0; i < n; i++) { // Find row w/max element in this int maxRow = i; // column, at or below diagonal for (int k = i + 1; k < n; k++) if (Math.abs(q[k][i]) > Math.abs(q[maxRow][i])) maxRow = k; if (maxRow != i) // If row not current row, swap for (int j = i; j < m; j++) { double t = q[i][j]; q[i][j]= q[maxRow][j]; q[maxRow][j]= t; q[maxRow][j]= t; } for (int j = i + 1; j < n; j++) { // Calculate pivot ratio double pivot = q[j][i] / q[i][i]; for (int k = i; k < m; k++) // Pivot operation itself q[j][k] -= q[i][k] * pivot; } } } 4Back Substitution private static void back_solve(double[][] q) { int n = q length; int n = q.length; int m= q[0].length; for (int p= n; p < m; p++) { // Loop over p columns for (int j = n - 1; j >= 0; j--) { // Start at last row double t = 0.0; // t- temporary for (int k = j + 1; k < n; k++) t += q[j][k] * q[k][p]; q[j][p]= (q[j][p] - t) / q[j][j]; } } } Variations Multiple right hand sides: augment Q, solve all eqns at once 33 112-2 55 77 87 87 2 4 3 35 75 -1 1 -3 0 -5 38 52 Matrix inversion (rarely done in practice) 313 1 -22 2 4 3 1 -3 0 100 1 0 0 0 1 0 0 0 1 ### # # # 0 # # 0 0 # @@@@ @@@ @ @@ @ @A I A-1 A • ? = I Q ? = A-1 56Invertpublic static double[][] invert(double[][] a) {int n = a.length; // Number of unknownsdouble[][] q = new double[n][n+n];for (int i = 0; i < n; i++) f (i t j 0 j j ) // F t ifor (int j = 0; j < n; j++) // Form q matrixq[i][j]= a[i][j];// Form identity matrix in right half of qfor (int i = 0; i < n; i++) q[i][n+i]= 1.0;forward_solve(q); // Do Gaussian eliminationback_solve(q); // Perform back substitutiondouble[][] x= new double[n][n]; // Extract R half of q for (int i = 0; i < n; i++) // which contains inverse for (int j= 0; j < n; j++)x[i][j]= q[i][j+n];return x;}// Method multiply() in download// Example use in GaussTest in downloadLU decomposition• We can write matrix A as the product of two matrices L and U:000luuuuaaaa• We can solve33323130222120111000000000llllllllll33232213121103020100000000uuuuuuuuuu.33323130232221201312111003020100aaaaaaaaaaaaaaaa=bxULxULxA=⋅⋅=⋅⋅=⋅)()(by first solving for a vector yand then solvingWhy? Solving each is trivial: forward, back substitutionbyL=⋅yxU=⋅Why and How • This is perhaps twice as fast as GaussianThis is perhaps twice as fast as Gaussian elimination (count steps) • L and U do not depend on b, so we can solve as many right hand sides as we wish • How: Crout’s method – We can decompose matrix A into matrices L and U by arranging the equations in a given orderarranging the equations in a given order – The rearrangement is subtle; we don’t cover it in class since you’ll never need to implement or modify it • Java implementation is on the Web site, based on Press et al, Numerical Recipes Class LU • Constructor: LU(double[][] a) – Stores LU decomposition in a single matrix Stores LU decomposition in a single matrix • All lii = 1.0 in matrix L • We store all u elements and all non-diagonal l elements in LU • Methods: – public double[] solve(double[] b) – public double[][] solve(double[][] b) – public double[][] inverse() – public double determinant()public double determinant() – public double[] improve(double[] b, double[] x) • See download for code and LUTest class for examples of usage – You can use it as a ‘black box’ – Use this in preference to class Gauss 78Other linear system algorithms• Banded matricesBanded matrices• Sparse matrices• Singular value decompositions (SVD)– Should be used in least squares computations• Cholesky decomposition (A= L LT)– Square, symmetric, positive definite matricesUdi ti– Used in econometrics• And others…– Almost all are based on pivot operationsLinear system model: Rail performance• Compute running time, including delays, for trains on a single track railroadon a single track railroad– Traffic in both directions: east and west – Three types of train (6 classes of train, including direction)Class D escription Velocity0WB way freight ‐251WB thru …