Microsoft Word - AERO 300 Laboratory 4 - Iterative Methods to Solve Matrix Equations.docx
AERO 300 Laboratory 4
Iterative Methods to Solve Linear Matrix Equations
The pre-lab assignment (Section 4) is due at the beginning of this lab (can be hand-written or
typed). The lab (Section 5) will be due before your next lab section (.zip file submitted to
polylearn).
1 Objectives
This lab introduces concepts and methods used to solving matrix equations iteratively. The
topics covered are:
• Solving linear systems of equations like !" = $ with MATLAB.
• Jacobi iterative method
• Gauss-Seidel iterative method
• Convergence criteria
At the completion of this lab you should be able to solve matrix equations using both the Jacobi
and Gauss-Seidel method.
2 Introduction
Linear systems of equations can be expressed in matrix form. Consider the system of equations:
3" + 5( = 3 (1)
2" + 6( = 1 (2)
This system can be expressed in the form !" = $ , where ! is a matrix of the coefficients, " is
the solution vector we seek, and $ is the right-hand side (or known part) of the system. Thus, the
above system is written in matrix form as:
XXXXXXXXXX. -
"(. = XXXXXXXXXX)
There are many methods and schemes for solving systems of equations including: alge
aic
substitution, alge
aic addition, graphically, matrix decomposition, Gaussian elimination, and
iterative schemes.
Iterative schemes differ from direct methods in that they converge on a numeric approximation
of the exact solution. Of course, direct and iterative schemes are always limited by machine
precision. So, normally an engineer should balance the scheme with the size of the problem.
For relatively small problems, direct methods are sufficient. When the size of the problem is
large (or the system is sparse – most of the elements of the matrix ! are zeros), iterative
schemes are prefe
ed. Practical uses of iterative schemes include applications in circuit
analysis, boundary value problems and partial differential equations. In Aerospace Engineering,
iterative schemes for solving systems of equations are used in Computational Fluid Dynamics
(CFD), Finite Element Analysis (FEA) for structures, and Computational Thermal Dynamics to
understand heat flow.
Iteration is a popular technique for finding roots of equations. A generalization of fixed-point
iteration can be applied to systems of linear equations to produce accurate results. The Jacobi
iteration method is attributed to Carl Jacobi XXXXXXXXXXand Gauss-Seidel iteration is
attributed to Johann Carl Friedrich Gauss XXXXXXXXXXand Philipp Ludwig von Seidel (1821-
XXXXXXXXXXIn both schemes, the solution vector " of the system of equations is found without
directly solving for the inverse of !.
2.1 Jacobi Iteration Method
Consider the following set of coupled linear equations. (Note that this method can be applied to
any number of equations. A set of four was chosen for this explanation.)
!" = $ (4)
Now, if we let ! = / XXXXXXXXXX, we can write
" = /234$ − (0 + 1)"8 (5)
Now, using an iterative scheme, we have,
"9:3 = /234$ − (0 + 1)"98 (6)
Jacobi iteration uses all previous step information to calculate the next approximate solution. In
the above equations, the left-hand side is the “next approximate solution” and the right-hand
side is the “previous step information.” Equation (6) is known as Jacobi Iteration and can used
to solve for the unknown vector " in equation (4). Since we do not have previous time step data
when we first begin the iteration, we must make an initial guess, ";. The better the guess, the
faster the iteration will converge.
As with some of our other methods, we are using MATLAB to “cheat” with built in matrix
multiplication. To implement the algorithm using the coefficients of ! and $ and for loops, we
can reformulate the problem by noticing, for example, that the first two elements of " are
updated as,
"39:3 = 1<33 4$3 − <3="=9 − <3>">9 −⋯− <3@"@98
"=9:3 = 1<== 4$= − <=3"39 − <=>">9 − ⋯− <=@"@98
(7)
This relationship can be written in general as,
"A9:3 = 1
@
DE3,AGD
H
(8)
2.2 Gauss-Seidel Method
Notice in equation (7) how the previous (or kth) value of "3 is used to find the new estimate of "= (i.e., at the k+1th step). Why not use the most recent estimate to calculate the estimate of
the solution? This idea is known as Gauss-Seidel Iteration. To formulate the Gauss-Seidel
Method, use the same factorization for the ! matrix, ! = / XXXXXXXXXXThen, the Gauss-Seidel
iteration scheme is given by
"9:3 = /234$ − 0"9:3 − 1"98
(9)
We can then solve for "9:3 as,
"9:3 = (0 + /)234$ − 1"98 (10)
Here again we are “cheating” by using MATLAB to take inverses and perform matrix
multiplication. We can use a for loop instead, if we re-write the IJK equation as,
"A9:3 = 1
A23
DE3
− C
@
DEA:3
H
(11)
2.3 A Note About Summations
Let us say we have the equation ( = 5 + ∑ 2MNDE3 . We could write this out in one line:
( = 5 + 2 ∙ 1 + 2 ∙ 2 + 2 ∙ 3 + 2 ∙ 4 + 2 ∙ 5
or we could add in each piece of the summation in steps:
( = 5
( = ( + 2 ∙ 1 = 7
( = ( + 2 ∙ 2 = 11
( = ( + 2 ∙ 3 = 17
( = ( + 2 ∙ 4 = 25
( = ( + 2 ∙ 5 = 35
giving us a final answer of 35 through either method. The summation works like a for loop, for
j=1:5, add 2*j to y. So the equation in MATLAB would be:
y = 5;
for j = 1:5
y = y + 2*j;
end
In general, summations can be implemented in MATLAB with for loops.
2.4 Notes on Matrix Multiplication
As this lab requires you to manipulate vectors and matrices, some notes and tips on indexing
matrices are provided here. Consider the matrix A:
! = R XXXXXXXXXX9U
In MATLAB, if we want to refer to the element in the first column and second row, we write:
A(2,1)
This will pull the value, 2, from matrix A. We can store this element in a variable called x by
writing:
x = A(2,1);
The variable x now has the value 2 stored in it. In general, if we want to refer to the information
stored in row i, and column j, we write:
A(i, j)
There are times when we will want to access an entire column of a matrix. This is done using the
colon operator. If we want to store the second column of matrix A into the vector y, we write:
y = A(:,2);
This statement will store every row in column two of the A matrix in y. If we now display y we
get:
disp(y)
5
6
6
Similarly, if we want to access the entire third row of A, we write:
A(3,:)
This statement retrieves every column in row three of matrix A. If only some subset of rows or
column is desired, then these can also be accessed using the colon operator. If we want the top
two values in column three, we would write:
x = A(1:2, 3);
Now, if we display x, we get:
disp(x)
1
8
Similarly, if we want the 2x2 matrix comprised of the lower right portion of matrix A we can write:
A(2:3,2:3)
ans =
XXXXXXXXXX
XXXXXXXXXX
In general, if we want to access rows, a to b, and columns c to d, we can write:
A(a:b, c:d)
3 Help
Potentially useful terms for this lab:
• linspace
• for
• function
• tic, toc
• while
• disp
• plot
• subplot
• axis
• grid
• hold
• legend
• title
• xlabel
• ylabel
4 Pre-lab Assignment
NOTE: This must be done before lab. If it is not completed, you will not be allowed into lab.
For a general 4 × 4 matrix ! = W
• write out all four equations for !" = $ similar to Equation (7).
• Starting from equation (9), show that that Gauss-Seidel method can be written out in
summation form as in equation (11)
5 Lab Assignment
Complete the following problems with a single script and custom functions as appropriate. Important lines
of code should include descriptive comments. You are encouraged to work in groups, however the code
you submit must be your own. All graphs/figures/sketches should include a grid, labels, legend (if
necessary), and the appropriate fontsize/linewidth/markersize.
1. Write pseudocode for Jacobi and Gauss-Seidel iterative methods using Equations (8) and (11).
Your algorithm should take A and b as inputs and return the solution x. Your algorithm should
also check for convergence before running the iterative scheme. If the convergence test does
not pass, warn the user. Continue with the iterative scheme, but stop after a fixed number of
iterations. If the solution does not converge, pass an e
or to the user and let them know the
solution failed.
Andy G Cristales
Andy G Cristales