Assignment 5 Due Date: March 11
A. Develop a routine power_iteration: int pow_iteration(int n,double **A,double shift, int max_it,double err,double *l,double *v); to find an eigenvalue l and eigenvector v of the n*n matrix A, by performing power iteration on the matrix B = A-shift*I to find the largest magnitude eigenvalue lmax of B, and then returning the corresponding eigenvalue l of A (i.e. lmax + shift) in *l. The routine should return 1 on success or 0 on failure (too many iterations, see below).
Here max_it is the maximum number of iterations to use and err is the required accuracy - iteration stops if successive iterates for the eigenvalue lmax differ by less than err. The eigenvalue of A (not B) is returned in *l and the eigenvector in v (v is an evector of both A and B). The vector v should be initialized to an initial guess (can be anything) before calling.
The code should look something like this pseudo-code where we seek lmax, the largest evalue of B.
B = A - shift*I; /* Allocate B and use copy_dmatrix() */ lmax_prev = HUGE; while (k < max_it) { k++; /* Compute w = B*v: */ w = B*v; /* Allocate w and Use apply_dmatrix(n,n,B,v,w); */
/* Find element of w with largest magnitude: */ for (i=1,lmax=w[0]; i<n; i++) if (fabs(w[i])>fabs(lmax)) lmax = w[i];
/* Normalize w and copy into v as the new iterate: */ for (i=0; i<n; i++) w[i] /= lmax; v = w; /* (use copy_dvector() */
/* Check for convergence */ if (fabs(lmax-lmax_prev) < err) break; lmax_prev = lmax; } if (k == max_it) return 0; *l = lmax+shift; /* We computed an eigenvalue of A-shift*I, so add shift to it */ return 1; }
You can verify the code by testing on the matrix A which is diagonal with elements 2 3 4 on the diagonal. Start v with initial guess (1 1 1) and set shift to 1.0. You should find 3 as the largest evalue of B and the routine should return 4 as the eigenvalue of A. Take max_it = 100 and err = 1.e-10
B. Develop a routine inverse_power_iteration: int inverse_pow_iteration(int n,double **A,double shift, int max_it,double err,double *l,double *v); to find an eigenvalue l and eigenvector v of the n*n matrix A, by performing power iteration on the matrix C = 1/(A-shift*I) to find the largest magnitude eigenvalue lmax of B, and then returning the corresponding evalue of A: l = 1/lmax + shift, in *l. The routine should return 1 on success or 0 on failure (too many iterations or no LU decomposition)
The code requires only three lines to be added to the previous code. Basically we need to replace the line w = B*v above by w = 1/B*v. This is equivalent to B*w = v. So instead of APPLYING B to v, we instead solve a linear equation with v as the RHS. To do this we get the LU decomp of B, store it (including the permutation p) and then at each iteration we use the LU to solve B*w = v.
The code should look something like this pseudo-code where we seek lmax, the largest evalue of 1/B.
/* Construct B=A-shift*I by copying A and substracting shift from its diagonal: */ B = A - shift*I /* Find and store the LU decomp of B using pivoting: */ if (find_lu_pp(n,B,p,LU,prt) == 0) return 0; /* LU decomp failed */
lmax_prev = HUGE; while (k < max_it) { k++; /* Compute w = 1/B **v : */ /* So equivalently solve B*w = v: */ solve_lu_pp(n,p,LU,v,w);
/* Find element of w with largest magnitude: */ for (i=1,lmax=w[0]; i<n; i++) if (fabs(w[i])>fabs(lmax)) lmax = w[i];
/* Normalize w: */ /* Normalize w and copy into v as the new iterate: */ for (i=0; i<n; i++) w[i] /= lmax; v = w;
if (fabs(lmax-lmax_prev) < eps) break; lmax_prev = lmax; } if (k == max_it) return 0; *l = shift + 1.0/lmax; /* since lmax was the evalue of 1/(A-shift) */ return 1;
You can verify the code by testing on the matrix A which is diagonal with elements 2 3 4 on the diagonal. Start v with initial guess (1 1 1) and set shift to 1.5. You should find 2.0 as the largest evalue of 1/B (i.e. .5 is the smallest evalue of B) and the routine should return 2 = 1.5 + 1/2.0 as the corresponding eigenvalue of A. Take max_it = 100 and err = 1.e-10.
C. Results
Let A be the nxn Hilbert matrix with elements A[i][j] = 1/(i+j+1), i,j>=0, with n = 8. Take max_it = 1000 and err = 1.e-14. Apply the routine in part A to find the largest eigenvalue of A by using shift = 0.0. Apply the routine in part B to find the smallest eigenvalue of A by using shift = 0.0 Apply the routine in part B to find the closest eigenvalue to .5 by using shift = .5
Input the three eigenvalues that you find on the web page for assignment results 5.
Note: you will find that the smallest evalue is VERY small. This is why the condition number of this matrix is so large. Condition number = ||A||*||1/A| = lmax(A)*lmax(1/A) = lmax(A)/lmin(A). This explains the bizarre huge errors we saw in Assignment 4 when solving equations with this matrix.