Let's first setup the matrices for Jacobi iteration. Type the following:
Make sure you understand each step here. diag(A) creates a vector whose elements are the diagonal elements of A. diag(diag(A)) then returns a matrix with those values on the diagonal. Notice A\b gave us the solution of Ax=b which we will try to reproduce in a moment.
The eig(LQ) command returned the eigenvalues of LQ. Notice some of them are not real.
Are they all in the unit circle in the complex plane?
What if they were?
Now type:
Here x is a matrix whose n-th column is the n-iterate x(n). Typing x' makes these interates appear in the n-th row (i.e., x' is the transpose of x). The plot command plots each column of x' and one can see the components of x(n) approach the solution.
Now let's do Gauss-Seidel iteration. The only difference is in how we define Q. Type:
Note how the "tril" extracts the lower triangular part of the matrix A.
Repeat the previous procedure to determine the first 15 iterates using Gauss-Seidel iteration using a seed value x(:,1)=[1;1;1] (same A and b). Your output should be:
To do the same for SOR you may need the "tril(A,n)" command. Read about it in the "helpwin". Try:
tril(A,-1)
tril(A,0)
tril(A,1)
to see what it does. Using this and the "diag" command you should be able to construct LQ and bQ for SOR Iteration. Use w=1.1 as the relaxation parameter and then perform the iterates with an initial seed value x(:,1)=[1;1;1] as before. Your output should be:
If we change the initial matrix A to (same b):
What do the LQ eigenvalues tell you about the convergence of SOR with w=1.1?
Test your hypothesis by performing the iteration using matlab.
Iteration using an M-file code and convergence rates:
Download the following M-file: Iterate.m
This code does Jacobi,Gauss-Seidel or SOR iteration depending on its inputs. Suppose you wanted to perform N=10 Jacobi iterations on the system
Ax = b
you'd type:
>> [x,ELQ]=Iterate(A,b,10,'jacobi',0);
>> x
The n-th row of the matrix x is the n-th iterate. Thus,
>> x(4,:)
is the 4rth iterate. Convergence (when it occurs) should appear as each column approaching some constant value. The command
>> plot(x)
would plot all the components of the n-th iterate as a function of n.
The vector ELQ are the eigenvalues of the matrix LQ. As in class, if these eigenvalues are all in the unit circle in the complex plane, the method will converge to the solution regardless of what initial guess is used (The code Iterate.m always uses the zero vector as an initial seed value).
The last argument in "Iterate" is the relaxation parameter for SOR. It must be assigned a value even for Jacobi and Gauss-Seidel iteration. Thus,
>> [x,ELQ]=Iterate(A,b,10,'jacobi',100000)
>> x
would have done the same thing as the commands above.
For 15 iterates using Gauss-Seidel, type:
>> [x,ELQ]=Iterate(A,b,15,'gauss',0)
>> x
For 17 iterates using SOR with a relaxation parameter w=1.23, type:
>> [x,ELQ]=Iterate(A,b,15,'sor',1.23)
>> x
If you look at the code for Iterate.m (using emacs) you will notice a new "if...then"-like programming structure. There the "switch....case....end" construct is used. You can read more about that under the helpwin like we did in previous computer labs.
Method Comparsion
With b=[1;2;3], execute the following:
Then type:
What you see are plots of (clockwise from top left figure) the components of iterate xn versus n in Jacobi, Gauss-Seidel and SOR, respectively. The graph should look like:
What feature of the matrix A should have made you suspect slow convergence in Jacobi iteration?
Power Method for eigenvalue of largest magnitude:
We need to know the matlab command for taking the norm of a vector:
>>b=[1;2;3];
>>bn=b/norm(b);
In the above, bn will be a unit vector in the direction of b since norm(b) is the norm of the vector.
Recall the Power Method (version discussed in class) had the following algorithm:
Try the associated matlab code for 10 iterations:
You will see the lambda values converge to 4.6432 in the screen output. Compare this with the eigenvalues of A computed using the "eig(A)" command. Then type:
>> A*y-lambda*y
What is the significance of all components of the resulting vector being small?
Now, download the following M-file: Power.m . Executing
>> [y,lambda]=Power(A,x1,N);
will do N Power method iterations for the matrix A using the initial guess x1. The vector lambda contains all the lambda values and y the associated eigenvector approximations. For example, type
>> x1=ones(3,1);
>> [y,lambda]=Power(A,x1,10);
>> plot(lambda)
will give a plot of the lambda values converging to the eigenvalue of A of largest absolute value.
Now, type:
>> A=[1 2 3;3 4 2;1 -1 2];
>> x1=ones(3,1);
>> [y,lam]=Power(A,x1,40);
>> lamtrue=eig(A);
>> error=abs(lam-lamtru(2));
>> plot(1:40,log(error),'*');
Here's a really hard question.
Why is the line straight and what is it's slope for a general matrix A?