Worksheet #3 - PHY301, Fall 2003

Iterative methods of linear algebra I

Due Friday Oct 31

The problems we solved last week are initial value problems in ordinary differential equations. The methods we used were explicit methods which calculate time updates using information at previous times. For boundary value problems and for partial differential equations, we often need more sophisticated methods called implicit methods. Before we can apply implicit methods however, we must learn to solve problems in linear algebra.

Two of the most common problems in linear algebra are finding solutions to a large set of linear equations and finding eigenvalues of large matrices. There are many standard computer packages to do these calculations (mathematica has some) however the majority are not optimized for physics problems. Physics problems usually involve sparse matrices which have most of their entries being zero. In that case it is very important to design algorithms which only operate with the non-zero elements of the matrix. Sparse linear algebra problems arise in solving partial differential equations using implicit methods as shall be illustrated in the next worksheet.

In this worksheet you will develop a C++ code for a simple iterative method for finding the largest eigenvalue of a real symmetric matrix. Next worksheet you will apply this method to find the lowest energy state of a quantum particle in a box.

Structuring of code

You should try to write code that is self-explanatory. So, the code you hand in for this assignment needs to have,

(i) Lots of comments that explain each step in the algorithm;
(ii) Variable names that are self-explanatory (e.g. do NOT use ``m'' to represent a matrix. You need to call it ``matrix'' or another clear descriptor).

Here is the ``iterative power method''

The eigenvalue of largest magnitude of a real symmetric square matrix can be found by taking an arbitrary vector and applying the matrix to it many times (this is sometimes called the ``iterative power method''). So if A is a matrix and u₀ is a (normalised) starting vector, then

$$A^n \vert u_0> \sim b \lambda_0^n \vert v_0>$$ (1)

where $\lambda_0$ is the eigenvalue of largest magnitude and v₀ is its associated eigenvector, and b is the projection <v₀|u₀>. It is quite simple to prove this using the expansion of |u₀> in the eigenvectors of A.
Problem

Write a c++ code for the iterative power method to find the eigenvalue of largest magnitude of a real, symmetric, square matrix (allow the size of the matrix to be a variable). Generate the matrix and the starting vector inside the program. You can choose the form of these. Note that the starting vector should not be ``special'' (e.g. don't use (1,0,0,0....) but (0.2343,0.1234,-0.453,...) is OK.

To find the largest eigenvalue iterate the following: Start with a normalized vector for |ulast>.
1. Multiply by the matrix: |uprod>=A|ulast>.
2. The approximate eigenvalue: eapprox = (<uprod|uprod>)^1/2. Print out eapprox - it is your approximation to the largest eigenvalue.
3. Generate the next |ulast>=|uprod>/eapprox. Return to 1. |ulast> is the approximate largest eigenvector.

Demonstrate that your algorithm works by checking it on a couple of matrices for which you know the largest eigenvalue. You can do this by finding the eigenvalues of your test matrices using Mathematica. You should also impose a convergence criterion on the eigenvalue, e.g. four digit accuracy. When writing your code, put the matrix product step and the inner product in two separate subroutines.

Hand in a hardcopy of your code, one of your test matrices and its eigenvalues (e.g. from Mathematica) and the output from your fortran code for this matrix (i.e. the sequence of estimates to eapprox). Plot the sequence of estimates on a log-linear plot. Note that at large iteration number the convergence is exponential. Why?