1.
00 Lecture 31
Systems of Linear Equations
Reading for next time: Numerical Recipes, pp. 129-139
http://www.nrbook.com/a/bookcpdf.php
Systems of Linear Equations
3x0 + x1 - 2x2 = 5
2x0 + 4x1 + 3x2 = 35
x0 - 3x1 = -5
3 1 -2 x0 5
2 4 3 x1 = 35
1 -3 0 x2 -5
A x = b
3 x 3 3 x 1 3 x 1
1
Algorithm to Solve Linear System
x0 b0
Create matrix =aij
x1 = b1
x2 b2
A x b
x0 b0 =aij
Forward solve 0 x1 = b1
0 0 x2 b2
A x b
x0 b0 x0
Back solve x1 b1
0 = x1
0 0 x2 b2 x2
A x b
Gaussian Elimination: Forward Solve
3 1 -2 5 Form Q for convenience
Do elementary row ops:
Q= 2 4 3 35 Multiply rows
1 -3 0 -5 Add/subtract rows
A b
Make column 0 have zeros below diagonal
Pivot= 2/3
Pivot= 1/3 3 1 -2 5
0 10/3 13/3 95/3 Row 1= row 1 - (2/3) row 0
0 -10/3 2/3 -20/3 Row 2= row 2 - (1/3) row 0
Make column 1 have zeros below diagonal
3 1 -2 5
Pivot= -1 0 10/3 13/3 95/3
0 0 15/3 75/3 Row 2= row 2 + 1 * row 1
2
Gaussian Elimination: Back Solve
3 1 -2 5
0 10/3 13/3 95/3
0 0 15/3 75/3 (15/3)x2=(75/3) x2 = 5
3 1 -2 5
0 10/3 13/3 95/3 (10/3)x1 + (13/3)*5= (95/3) x1 = 3
0 0 15/3 75/3
3 1 -2 5 3x0 + 1*3 - 2*5 = 5 x0 = 4
0 10/3 13/3 95/3
0 0 15/3 75/3
A Complication
0 1 -2 5
2 4 3 35 Row 1= row 1 - (2/0) row 0
1 -3 0 -5
Exchange rows: put largest pivot element in row:
2 4 3 35
0 1 -2 5
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.
3
Gaussian Elimination
// In class Matrix, add:
public static Matrix gaussian(Matrix a, Matrix b) {
int n = a.data.length; // Number of unknowns
Matrix q = new Matrix(n, n + 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) // Form q matrix
q.data[i][j]= a.data[i][j];
q.data[i][n]= b.data[i][0];
}
forward_solve(q); // Do Gaussian elimination
back_solve(q); // Perform back substitution
Matrix x= new Matrix(n, 1);
for (int i = 0; i < n; i++)
x.data[i][0]= q.data[i][n];
return x;
}
Forward Solve
private static void forward_solve(Matrix q) {
int n = q.data.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.data[k][i]) > Math.abs(q.data[maxRow][i]))
maxRow = k;
if (maxRow != i) // If row not current row, swap
for (int j = i; j <= n; j++) {
double t = q.data[i][j];
q.data[i][j]= q.data[maxRow][j];
q.data[maxRow][j]= t;
}
for (int j = i + 1; j < n; j++) { // Calculate pivot ratio
double pivot = q.data[j][i] / q.data[i][i];
for (int k = i; k <= n; k++) // Pivot operation itself
q.data[j][k] -= q.data[i][k] * pivot;
}
}
}
4
Back Solve
private static void back_solve(Matrix q) {
int n = q.data.length;
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.data[j][k] * q.data[k][n];
q.data[j][n]= (q.data[j][n] - t) / q.data[j][j];
}
}
// GaussTest is in your download for simple tests of this code
Variations
Multiple right hand sides: augment Q, solve all eqns at once
3 1 -2 5 7 87
2 4 3 35 75 -1
1 -3 0 -5 38 52
Matrix inversion:
3 1 -2 1 0 0 # # # @ @ @
2 4 3 0 1 0 0 # # @ @ @
1 -3 0 0 0 1 0 0 # @ @ @
A I A-1
A•?=I
Q
? = A-1
5
Exercise
• Download GElim and Matrix
• Compile and run GElim:
Replace row
with result
© Oracle. All rights reserved. This content is excluded from our Creative
Commons license. For more information, see http://ocw.mit.edu/fairuse.
Exercise
• Experiment with the following 3 systems:
– Use pivot, back subst, divide on selection, etc. not solve
System 1: The 3x3 matrix example in the previous slides. Click
on Demo to load it.
System 2:
4 6 -3 10
2 5 9 12
8 8 -27 6
System 3:
12 -13.5 3 0.5 2.75
8 -9 4 2.5 3.5
3 6 1.5 2 4.25
2 1.5 4 12 6
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
Heat Transfer Exercise
60
x0 x1 x2 x3 4 by 4 grid of points
x x5 x6 x7 on the plate produces
80 4 40
x8 x9 x10 x11 16 unknown temperatures
x12 x13 x14 x15 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
7
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
Ax= b
j
A x b
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 x0 140
1 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 x1 60
2 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 x2 60
3 0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0 x3 100
4 -1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0 x4 80
5 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 x5 0
6 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 x6 0
i 7 0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0 x7 40
8 0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0 x8 80
9 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 x9 0
10 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 x10 0
11 0 0 0 0 0 0 0 -1 0 0 -1 4 0 0 0 -1 x11 40
12 0 0 0 0 0 0 0 0 -1 0 0 0 4 -1 0 0 x12 100
13 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 x13 20
14 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 x14 20
15 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 x15 60
8
Heat Transfer System
16
a00 a01 a02 a0,15 x0 b0
a10 a11 a12 a1,15 x1 b1
16 a20 a21 a22 a2,15 x2 = b2
a15,0 a15,1 a15,15 x15 b15
A x b
Contains 0, -1, 4 Known temperatures
coefficients in (often 0 but use edge
(simple) pattern temperatures when close)
16 unknown interior temperatures
Heat Transfer Result
60
66 59 55 50
66 56 50 45
80 40
62 50 44 41
50 38 34 34
20
9
Exercise: Ax= b: Fill in green, orange
i=j j=i+1 j= ___
_______ j
A x b
j=i-1 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 x0 140
1 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 x1 60 N
j=___ 2 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 0 0 0 x2 60
3 0 0 -1 4 0 0 0 -1 0 0 0 0 0 0 0 0 x3 100
____ 4 -1 0 0 0 4 -1 0 0 -1 0 0 0 0 0 0 0 x4 80
5 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 0 x5 0
____ 6 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 0 0 0 x6 0 W
i 7 0 0 0 -1 0 0 -1 4 0 0 0 -1 0 0 0 0 x7 40
8 0 0 0 0 -1 0 0 0 4 -1 0 0 -1 0 0 0 x8 80
9 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 0 x9 0 E
10 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 0 -1 0 x10 0
11 0 0 0 0 0 0 0 -1 0 0 -1 4 0 0 0 -1 x11 40
12 0 0 0 0 0 0 0 0 -1 0 0 0 4 -1 0 0 x12 100
13 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 0 x13 20 S
14 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 -1 x14 20
15 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 4 x15 60
Heat Transfer Exercise, p.1
public class Heat { // Problem generator
public static void main(String[] args) {
double Te= 40.0, Tn=60.0, Tw=80.0, Ts=20.0; // Edge temps
int col= 4; // Te – east, Tn – north, Tw- west, Ts-south
int row= 4;
int n= col * row;
Matrix a= new Matrix(n,n);
for (int i=0; i < n; i++)
for (int j=0; j < n; j++) {
if (i==j) // Diagonal element (yellow)
a.setElement(i, j, 4.0);
else if (…)
// Complete this code in step 1:
// Green elements (4, or col, away from diagonal)
// Blue elements (1 away from diagonal)
// Set blue and skip orange where we go to
// the next row on the actual plate
// Relate i and j to determine blue, orange cells
// using the diagram on the previous slide
// Continued on next slide
10
Heat Transfer Exercise, p.2
Matrix b= new Matrix(n, 1); // Known temps
for (int i=0; i < n; i++) {
if (i < col) // Next to north edge
b.incrElement(i, 0, Tn); // incrElement, not setElement
if (…) // Step 2
// Complete this code for the other edges; no elses
// Add edge temperature to b; you may add more than one
// Look at the Ax=b example slide to find the pattern
// Use i, col, row to determine cells at the edge
}
Matrix x= Matrix.gaussian(a, b); // Problem solution
System.out.println("Temperature grid:"); // Output generator
for (int i=0; i< row; i++) {
for (int j=0; j < col; j++)
System.out.print(Math.round(x.getElement((i*row+j),0)+ ");
System.out.println();
}
}
}
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 (variable) 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. Youll
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, youll
find its not a solution. (Right sides dont quite equal left sides.)
11
MIT OpenCourseWare
http://ocw.mit.edu
1.00 / 1.001 / 1.002 Introduction to Computers and Engineering Problem Solving
Spring 2012
For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.
