Worksheet #3 - PHY301, Fall 2005

Iterative methods of linear algebra I

Due Friday Oct 28, 2005 - three weeks

The problems you solved in the last set were initial-value problems in ordinary differential equations. The methods used were explicit methods which calculate time updates using information from previous times. For boundary value problems and for partial differential equations, we often need more sophisticated methods called implicit methods. But before we can apply implicit methods, 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 often involve sparse matrices which have most of their entries being zero. In that case it is important to design algorithms that only operate with the non-zero elements of the matrix. Sparse linear algebra problems arise in solving partial differential equations using implicit methods as will 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. In the 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 strive to write code that is self-explanatory. So, the code you hand in for this assignment needs to have,
(i) Comments that explain each step in the algorithm;
(ii) Variable names that are descriptive -- e.g., do NOT just use ``m'' to represent a matrix.

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 nearly arbitrary 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. You can specify the dimension of the matrix inside the program by something like
    const int ndim=10;
    double Matrix[ndim][ndim];
or you can specify a maximum dimension by something like
    const int nmax=100;
    double Matrix[nmax][nmax];
and let your "user" specify the actual dimension to be used, with ndim <= nmax.

To find the largest eigenvalue, use the following steps:
  1. Create an "arbitrary" vector |ulast>.
Note that this starting vector should not be ``special'' --- e.g. don't use (1,0,0,0....) --- because you could be unlucky enough to find that it has no component along the desired eigenvector. A bunch of randomly made-up values like (0.2343, 0.1234, -0.4531,...) will be safe.
  2. Normalize that vector by dividing each component by the square root of <ulast|ulast>
After you do that, you will have <ulast|ulast>=1.
  3. Multiply by the matrix: |uprod>=A|ulast>.
  4. Calculate the approximate eigenvalue: eapprox = <uprod|ulast>. Print out eapprox.
  5. Generate the next (unnormalized) eigenvector approximation |ulast>=|uprod>.
  6. Quit if the iteration has converged. (A nice way to tell when the iteration has converged is take the dot product of the normalized ulast> with the previous normalized ulast> --- the iteration has converged when this quantity is within, say, 10^(-12) of 1 or -1.
  7. Return to step 2.
On exit, |ulast> is the approximate eigenvector corresponding to the eigenvector with the largest absolute value.

When writing your code, create two separate functions to handle the matrix product and the inner product.

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.

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?