lecture_23_s2005

lecture_23_s2005 - 1.00 Lecture 23 Systems of Linear...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: 1.00 Lecture 23 Systems of Linear Equations Reading for next time: Numerical Recipes, pp. 129-139 Linear Systems a00x0 + a01x1 + a02x2 + … + a0,n-1xn-1 = b0 a10x0 + a11x1 + a12x2 + … + a1,n-1xn-1 = b1 … am-1,0x0 + am-1,1x1 + am-1,2x2 + … + am-1,n-1xn-1 = bm-1 • If n=m, we try to solve for a unique set of x. Obstacles: – If any row (equation) or column (variables) is linear combination of others, matrix is degenerate or not of full rank. No solution. Your underlying model is probably wrong; you’ll need to fix it. – If rows or columns are nearly linear combinations, roundoff errors can make them linearly dependent during computations. We’ll fail to find a solution, even though one may exist. – Roundoff errors can accumulate rapidly. While you may get a solution, when you substitute it into your equation system, you’ll find it’s not a solution. (Right sides don’t quite equal left sides.) • Large linear systems tend to be close to singular (degenerate). Beware! 1 Systems of Linear Equations 3x0 + x1 - 2x2 = 5 2x0 + 4x1 + 3x2 = 35 = -5 x0 - 3x1 31 24 1 ­3 A 3x ­2 3 0 3 x0 5 x1 = 35 ­5 x2 x=b 3x1 3x1 Algorithm to Solve Linear System x0 x1 x2 Create matrix = b0 b1 b2 A Back solve 0 0 0 0 0 x0 x1 x2 A’ Forward solve x x b’ 0 x0 x1 x2 b’0 b’1 b’2 A’ x =aij b = = b’0 b’1 b’2 =a’ij x0 x1 x2 b’ 2 Gaussian Elimination: Forward Solve 3 2 1 Q= 1 4 -3 -2 3 0 5 35 -5 b A Form Q for convenience Do elementary row ops: Multiply rows Add/subtract rows Make column 0 have zeros below diagonal Pivot= 2/3 Pivot= 1/3 3 0 0 1 10/3 -10/3 -2 13/3 2/3 5 95/3 -20/3 Row 1’= row 1 - (2/3) row 0 Row 2’= row 2 - (1/3) row 0 Make column 1 have zeros below diagonal Pivot= 1 3 0 0 1 10/3 0 -2 13/3 15/3 5 95/3 75/3 Row 2’’= row 2’ + 1 * row 1 Gaussian Elimination: Back Solve 3 0 0 -2 1 10/3 13/3 15/3 0 5 95/3 75/3 3 0 0 -2 1 10/3 13/3 15/3 0 5 95/3 75/3 3 0 0 -2 1 10/3 13/3 15/3 0 5 95/3 75/3 (15/3)x2=(75/3) x2 = 5 (10/3)x1 + (13/3)*5= (95/3) x1 = 3 3x0 + 1*3 - 2*5 = 5 x0 = 4 3 A Complication 0 2 1 1 4 -3 -2 3 0 5 35 -5 Row 1’= row 1 - (2/0) row 0 Exchange rows: put largest pivot element in row: 2 0 1 4 1 -3 3 -2 0 35 5 -5 Do this as we process each column. If there is no nonzero element in a column, matrix is not full rank. Gaussian Elimination // In class Matrix, add: // In class Matrix, add: static Matrix Matrix b) public static Matrix gaussian(Matrix a, Matrix b) { a.data.length; Number of unknowns int n = a.data.length; // Number of unknowns new Matrix(n, 1); Matrix q = new Matrix(n, n + 1); (int 0; n; i++) for (int i = 0; i < n; i++) { (int 0; n; j++) for (int j = 0; j < n; j++) q.data[i][j]= a.data[i][j]; q.data[i][j]= a.data[i][j]; q.data[i][n]= b.data[i][0]; q.data[i][n]= b.data[i][0]; } forward_solve(q); forward_solve(q); back_solve(q); back_solve(q); Form matrix // Form q matrix Do Gaussian elimination // Do Gaussian elimination Perform back substitution // Perform back substitution x= new Matrix(n, 1); Matrix x= new Matrix(n, 1); for (int i = 0; i < n; i++) (int 0; n; i++) q.data[i][n]; x.data[i][0]= q.data[i][n]; x; return x; } 4 Forward Solve private static void private static void forward_solve(Matrix q) { int n = q.data.length; q.data.length; for (int i = 0; i < n; i++) { // Find row w/max element in this (int 0; n; i++) Find row w/max element in this i; column, at or below diagonal int maxr = i; // column, at or below diagonal (int 1; n; k++) for (int k = i + 1; k < n; k++) (Math.abs(q.data[k][i Math.abs(q.data[k][i]) Math.abs(q.data[maxr][i])) if (Math.abs(q.data[k][i]) > Math.abs(q.data[maxr][i])) k; maxr = k; (maxr i) If row not current row, swap if (maxr != i) // If row not current row, swap (int i; <= n; j++) for (int j = i; j <= n; j++) { q.data[i][j]; double t = q.data[i][j]; q.data[i][j]= q.data[maxr][j]; q.data[i][j]= q.data[maxr][j]; q.data[maxr][j]= t; q.data[maxr][j]= t; } for (int j = i + 1; j < n; j++) { // Calculate pivot ratio[i] ata[i][i]; double pivot = q.data[j][i] / q.data[i][i]; for (int k = i; k <= n; k++) // Pivot operation itself [k] q.data[i][k] q.data[j][k] ­= q.data[i][k] * pivot; } } } Back Substitution private static void private static void back_solve(Matrix q) { int n = q.data.length; q.data.length; for (in j = n ­ 1; j >= 0; j­­) { (int j­­) // Start at last row ­­ 0.0; t­ double t = 0.0; // t­ temporary (int 1; n; k++) for (int k = j + 1; k < n; k++) += q.data[j][k] q.data[k][n]; t += q.data[j][k] * q.data[k][n]; q.data[j][n]= (q.data[j][n q.data[j][n] q.data[j][j]; q.data[j][n]= (q.data[j][n] ­ t) / q.data[j][j]; } } 5 Test Program import javax.swing.*; import javax.swing.*; class public class GaussTest { static void main(String args) public static void main(String args) { j; int i, j; term; double term; input JOptionPane.showInputDialog( String input = JOptionPane.showInputDialog( "Enter number of unknowns"); number of unknowns"); Integer.parseInt(input); int n = Integer.parseInt(input); new Matrix(n, n); Matrix a = new Matrix(n, n); new Matrix(n, 1); Matrix b = new Matrix(n, 1); Enter matrix // Enter matrix a (i 0; a.getNumRows(); i++) for (i = 0; i < a.getNumRows(); i++) (j 0; a.getNumCols(); j++) for (j = 0; j < a.getNumCols(); j++) { =JOptionPane.showInputDialog JOptionPane.showInputDialog( input =JOptionPane.showInputDialog( a[" "][" "]"); "Enter a[" + i + "][" + j + "]"); Double.parseDouble(input); term = Double.parseDouble(input); a.setElement(i, j, term); a.setElement(i, j, term); Continued on next slide } // Continued on next slide Test Program, p.2 // Enter vector as 1­ matrix // Enter vector b as 1­column matrix (i 0; b.getNumRows(); i++) for (i = 0; i < b.getNumRows(); i++) { JOptionPane.showInputDialog( input = JOptionPane.showInputDialog( b[" "]"); "Enter b[" + i + "]"); Double.parseDouble(input); term = Double.parseDouble(input); b.setElement(i, 0, term); b.setElement(i, 0, term); } x= Matrix.gaussian(a, b); Solve it! Matrix x= Matrix.gaussian(a, b); // Solve it! System.out.println("Matrix a:"); a.print(); a.print(); System.out.println("Vector b:"); b.print(); b.print(); System.out.println("Solution vector x:"); x:"); x.print(); x.print(); } } 6 Variations Multiple right hand sides: augment Q, solve all eqns at once 3 2 1 -2 3 0 1 4 -3 5 35 -5 7 75 38 87 -1 52 Matrix inversion (rarely done in practice) 3 1 -2 243 1 -3 0 A 1 0 0 0 1 0 0 0 1 # 0 0 # # 0 # # # @ @ @ @ @ @ @ @ @ A-1 I Ax=b x= A-1 b Q Exercise • Download GElim and Matrix • Compile and run GElim: 7 Exercise • Experiment with the following 3 systems: – Use pivot, back subst, divide on selection, etc. not solve – The 3x3 matrix example in the previous slides. Click on “Demo” to load it. 4 2 8 12 8 3 2 6 5 8 -3 9 -27 -13.5 -9 6 1.5 10 12 6 3 4 1.5 4 0.5 2.5 2 12 2.75 3.5 4.25 6 Using Linear Systems • A common pattern in engineering, scientific and other analytical software: – Problem generator (model, assemble matrix) • Customized to specific application (e.g. heat transfer) • Use matrix multiplication, addition, etc. – Problem solution (system of simultaneous linear equations) • Usually “canned”: either from library or written by you for a library – Output generator (present result in understandable format) • Customized to specific application (often with graphics, etc.) • We did a pattern earlier: model-view-controller 8 Heat Transfer Example 60 x0 x 80 4 x8 x12 x1 x5 x9 x13 x2 x3 x6 x7 40 x10 x11 x14 x15 4 by 4 grid of points on the plate produces 16 unknown temperatures x0 through x15 20 T= (Tleft + Tright + Tup + Tdown )/4 Edge temperatures are known; interior temperatures are unknown This produces a 16 by 16 matrix of linear equations Heat Transfer Equations • Node 0: x0= (80 + x1 + 60 + x4)/4 4x0- x1 –x4= 140 • Node 6: x6= (x5 + x7 + x2 + x10)/4 4x6 –x5 –x7 –x2 –x10= 0 • Interior node: xi= (xi-1 + xi+1 + xi-n + xi+n )/4 4xi – xi-1 –xi+1 –xi-n –xi+n= 0 Node 0 1 2 3 4 5 6 7 0 4 -1 0 0 -1 0 0 0 1 -1 4 -1 0 0 -1 0 0 2 0 -1 4 -1 0 0 -1 0 3 0 0 -1 4 0 0 0 -1 4 -1 0 0 0 4 -1 0 0 5 0 -1 0 0 -1 4 -1 0 6 0 0 -1 0 0 -1 4 -1 7 0 0 0 -1 0 0 -1 4 9 Ax= b j i 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 1 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 2 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 3 0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0 4 -1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0 5 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 6 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 7 0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0 8 0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0 9 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 10 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 11 12 13 14 15 00000 00000 00000 00000 00000 00000 00000 -1 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 -1 0 0 -1 0 4 0 0 0 -1 0 4 -1 0 0 0 -1 4 -1 0 0 0 -1 4 -1 -1 0 0 -1 4 x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 140 60 60 100 80 0 0 40 80 0 0 40 100 20 20 60 Heat Transfer System 16 16 a00 a01 a02 … a0,15 a10 a11 a12 … a1,15 a20 a21 a22 … a2,15 … a24,0 a24,1 … a24,15 x0 x1 x2 … x15 A x Contains 0, -1, 4 coefficients in (simple) pattern b0 b1 = b2 … b15 b Known temperatures (often 0 but use edge temperatures when close) 16 unknown interior temperatures 10 Heat Transfer Result 60 66 66 80 62 50 59 56 50 38 55 50 44 34 50 45 41 34 40 20 Heat Transfer Exercise, p.1 public class Heat Problem generator public class Heat { // Problem generator static void main(String args) public static void main(String args) { Te= 40.0, Tn=60.0, Tw=80.0, Ts=20.0; Edge temps double Te= 40.0, Tn=60.0, Tw=80.0, Ts=20.0; // Edge temps col= 4; int col= 4; 4; int row= 4; row; int n= col * row; a= new Matrix(n,n); Matrix a= new Matrix(n,n); (int n; i++) for (int i=0; i < n; i++) (int n; j++) for (int j=0; j < n; j++) { (i==j) Diagonal element (yellow) if (i==j) // Diagonal element (yellow) a.setElement(i, j, 4.0); a.setElement(i, j, 4.0); else if (…) if (…) // Complete this code: Complete this code: // Green elements (4, or ncols, away from diagonal) Green elements (4, or ncols, away from diagonal) // Blue elements (1 away from diagonal) Blue elements (1 away from diagonal) Set blue and skip orange where we go to point // Set blue and skip orange where we go to a point on the next row on the actual plate // on the next row on the actual plate Continued on next slide // Continued on next slide 11 Heat Transfer Exercise, p.2 Matrix b= new Matrix(n, 1); Known temps Matrix b= new Matrix(n, 1); // Known temps n; i++) for (int i=0; i < n; i++) { (int (i col) Next to north edge if (i < col) // Next to north edge b.setElement(i, 0, b.getElement(i,0)+Tn); b.setElement(i, 0, b.getElement(i,0)+Tn); (…) if (…) Complete this code for the other edges; no ‘elses elses’! // Complete this code for the other edges; no ‘elses’! Add edge temperature to b; you may add more than one // Add edge temperature to b; you may add more than one Look at the Ax=b example slide to find the pattern // Look at the Ax=b example slide to find the pattern } Problem solution x= Matrix.gaussian(a, b); Matrix x= Matrix.gaussian(a, b); // Problem solution Output generator System.out.println("Temperature grid:"); // Output generator i< row; i++) for (int i=0; i< row; i++) { (int (int col; j++) for (int j=0; j < col; j++) System.out.print(Math.round(x.getElement((i*row+j),0)+“ "); System.out.print(Math.round(x.getElement((i*row+j),0)+“ "); System.out.println(); System.out.println(); } } } Solution with 10 by 10 grid 60 80 69 69 72 74 74 73 72 70 67 62 50 65 67 68 68 67 65 62 57 50 38 62 63 63 63 61 59 55 50 43 33 60 60 60 59 57 54 50 45 38 30 59 58 57 55 53 50 46 41 35 28 58 56 55 52 50 47 43 39 33 27 57 55 52 50 48 45 41 37 32 26 56 53 50 48 45 43 40 37 32 26 54 50 47 45 44 42 40 37 33 28 50 46 44 43 42 41 40 38 35 31 40 20 Same code as exercise; just change row= col= 10 You could create a grid of colors as Swing output If you use 100 by 100 grid, you’ll get very nice results 12 Other Applications • Solve systems with 1,000s or millions of linear equations or inequalities – Networks, mechanics, fluids, materials, economics – Often linearize systems in a defined range b a – Routines in this lecture are ok for a few hundred equations • They aren’t very good at picking up collinear systems. Check first-see Numerical Recipes – Otherwise, see Numerical Recipes 13 ...
View Full Document

This note was uploaded on 11/29/2011 for the course CIVIL 1.00 taught by Professor Georgekocur during the Spring '05 term at MIT.

Ask a homework question - tutors are online