CS 314 Notes, Week 3
Gaussian Elimination in Matlab
%The following is code for an m-file to do Gaussian elimination. To call the function, type "gauss(A,B)" where A and B are defined in the equation Ax=b
function x=gauss(A,B)
%A is NxN, B is Nx1. Solution X will be NX1
%Get dimensions of A
[N,N]=size(A)
%Initialize matric X with 0s
X=zeros(N,1)
%Create augmented matrix
Aug=[A B]
%In Gaussian elimination we divide the row of the augmented matrix, starting at the pivot, by the pivot
for p=1:N
piv=Aug(p,p)
for(k=p:N+1)
Aug(p,k)=Aug(p,k)/piv
end
%Note: shorthand notation would have been
%Aug(p,p:N+1)=Aug(p,p:N+1)/piv
%Next we make zero the elements below the pivot in the pivotal column
for (k=p+1:N) %for all the remaining rows
m=Aug(k,p)
for (i=P:N+1)
Aug(k,i)=Aug(k,i)-m*Aug(p,i)
end
end
end
%Note that this method does not work if the pivot is 0; we would have to swap rows.
%Now we have an upper triangular matrix:
%Now can use back substitution to solve:
for (k=N:1)
%To solve the example matrix, we'd set x3=c34/c33,x2=(c24-x3*x23)/c22, and x1=(c14-c13*x3-c12*x2)/c11
sum=0
for (i=k+1:N)
sum=sum+Aug(k,i)*x(i,1)
end
x(k,1)=(Aug(k,N+1)-sum)/Aug(k,k)
end
%So if call gauss(A,B) result will be X
%To reduce computing error you may want to switch rows so you are working with the row that has the largest pivots of the remaining rows:
An example using Gaussian elimination:
Solve the system
5x+3y+2z=0
2x+4y+7z=13
3x+2y+6z=11
Write the augmented matrix:
Divide 1st row by the pivot, 5:
Row 2=Row2-2*Row1:
Row3=Row3-3*Row1:
Divide Row 2 by pivot:
Row3=Row3-0.2*Row2
Now can solve using backwards substitution.
z=1
y=(3.1243-2.2143)/1=1
x=(2-0.4*1-0.6*1)/1=1
LU Factorization
It is a method used when you have multiple systems of equations:
Ax1=b1 Ax2=b2 Ax3=b3......
A has a triangular factorization if it can be expressed as L*U, where L is a lower triangular matrix and U is an upper triangular matrix:
We can rewrite the equation
Ax=B
in terms of the matrices L and U:
(LU)x=B
Then define Y=UX => LY=B
Once you obtain matrices L and U, there are two steps.
Step 1: Solve for Y in LY=B using forward substitution.
Step 2: Solve for x in Y=Ux using backward substitution.
L and U are obtained using a form of Gaussian elimination.
Suppose we have the following systems of equations:
6x1+x2-4x3=3
5x1+3x2+2x3=21
x1-4x2+3x3=10
6z1+z1-4z1=3
5z1+3z2+2z3=10
z1-4z2+3z3=10
Note that the A matrix is the same for both systems.
It follows from properties of the identity matrix that
We want the identity matrix to become lower triangular, and the matrix to its right to become upper triangular.
To do this, use Gaussian elimination on right matrix and place multiplier used in Gaussian Elimination in left matrix, in the same column as the pivot and in the row you are currently changing. For the example system, we do the following:
Skip the part where we make pivot 1, and set Row2=Row2-5/6*Row1. Then place m=5/6 in the proper place.
Row3=Row3-1/6*Row1, and m=1/6:
Row3=Row3-(-1.9231*Row2), and m=-1.9231:
Now the left matrix is lower triangular and the right matrix is upper triangular.
Note that LU must equal A, or something is wrong!
Step 1: Solve LY=Ux
y1=3, y2=18.5, y3=(10-1/6*3-(-1.9231*18.5))/1=45.07
Step 2: Solve Ux=Y
x3=45.0769/(13.9231)=3.2376
x2=(18.5-5.3333*3.2376)/(2.1667)=.5695
x3=(3-(-4)*3.2376-1*.5695)/6=2.5635
We can use the same L and U to solve the second system of equations, or any system that shares the same A.
Note that if pivots become 0, we must switch rows and add a middle matrix R that tells us the rows were switched.
When we have many systems with the same A, LU factorization is much more efficient than solving the systems separately:
Gaussian elimination is O(N^3) (because we have three nested loops)
Back substitution and forward substitution are each O(N^2) (because ((N-1)+(N-2)+....+1)=(N-1)(N-2)/2=O(N^2))
Adding these together, the efficiency is O(N^3).
With LU factorization we cut out the Gaussian elimination part.
To quantify the difference, suppose we have M systems of equations with the same A, where A can be factored into LU.
Using LU factorization, the problem is O(N^3)+M*O(N^2), whereas if we solved independently it would be M*O(N^3).
Iterative Methods for Systems of Linear Equations
These methods are faster than Gaussian elimination for sparse matrices (formally, matrices where O(N) elements are non-zero).
Jacobi Iteration
Suppose we have the following system:
1) 3*x+y=5
2) x+3*y=7
We can write:
3) x(k+1)=(5-y(k))/3
4) y(k+1)=(7-xk)/3
x(0)=0,y(0)=0
| k | xk | yk |
| 0 | 0 | 0 |
| 1 | (5-0)/3=5/3 | (7-0)/3=7/3 |
| 2 | (5-7/3)/3=.889 | (7-5/3)/3=1.77 |
| 3 | 1.074 | 2.037 |
| 4 | .987 | 1.975 |
This is often useful when finding numerical solutions of partial differential equations, when there can be as many as 100,000 variables.
Gaussian elimination would be O(100,000^3)=O(10^15)
Jacobi iterations would be much more eficient. Assume 1000 iterations, and that there are 3 non-zero coefficients in each row.
1000*3*100,000=3*10^8 -- much less than 10^15.
An even faster method is the Gauss-Seidel method.
Note that in Jacobi iterations, we were setting y(k)=(7-x(k))/3 even though we had already calculated x(k+1).
The Gauss Seidel equations for this system would be:
x(k+1)=(5-y(k))/3
y(k+1)=(7-x(k+1))/3
The difference is we are using values in the iteration as they are obtained.
| k | xk | yk |
| 0 | 0 | 0 |
| 1 | (5-0)/3=5/3=1.667 | (7-1.667)/3=1.77 |
| 2 | (5-1.77)/3=.1.074 | (7-1.074)/3=1.975 |
| 3 | (5-1.975)/3=1.008 | (7-1.008)/3=1.99 |
The convergence is faster.
The Jacobi and Gauss-Seidel methods do not always converge!
Convergence can be assured if A is strictly diagonal dominant.
A matrix A of dimensions NxN is strictly diagonal dominant if
ie, the absolute value of the element in the diagonal is greater than the sum of the absolute value of all the other elements in that row.
For example, in the following matrix,
for row 1 is must be true that
|a11|>|a12|+|a13|+|a14|+....+|a1N|
It is guaranteed that the Gauss-Seidel and Jacobi methods will converge is A is strictly diagonal dominant.
The system
3x+y=5
x+3y=7
is SDD so it converges.
The system
x+3y=5
3x+y=7
is not SDD so it is not guaranteed to converge. Applying the Jacobi method, we see
x(k+1)=5-3*y(k)
y(k+1)=7-3*x(k)
| k | xk | yk |
| 0 | 0 | 0 |
| 1 | 5-3*0=5 | 7-3*0=7 |
| 2 | 5-3*7=-16 | 7-3*5=-8 |
| 3 | 5-3*-8=29 | 7-3*-16=55 |
| 4 | 5-3*55=-160 | 7-3*29=-80 |
| 5 | 5-3*-80=245 | 7-3*-160=487 |
which is diverging!
The Gauss-Seidel method would diverge faster.
Newton's Method for Systems of Equations
Suppose we have
f1(x,y)=x^2-y-0.2=0
f2(x,y)=y^2-x-.3=0
To solve, we can extend the Newton method to equations of two variables.
Assume the Taylor expansion for two variables around x0 and y0:
Or, in matrix form:
(1)
The 2x2 matrix with partial derivatives is called the Jacobian.
We can use this equation to solve for [x1,y1]:
The inverse Jacobian is difficult to compute. Better to solve (1) for x(k+1),y(k+1).
Returning to the original example,
f1(x,y)=x^2-y-0.2=0
f2(x,y)=y^2-x-.3=0
we solve for x(k+1) and y(k+1):
Let's check the accuracy of x2 and y2 as solutions to the sysolution.
f1(x2,y2)=-.2947^2+.1132-.2=4.8e-5
f2(x2,y2)=-.1132^2-(-.2947)-.3=0.007
Pretty close to (0,0).