Description
Your task. Set δ = 10−8
. Run the experiment described in example Ex5 first. Then, for each of the
four matrices A from examples Ex1, Ex2, Ex3, and Ex4, take n = 42, and compute the eigenvalues
of A + δQ for 100 random orthogonal matrices. Plot the eigenvalues of all 100 computations
on one plot. If the vector zeig contains all your values, you might want to plot them with the
command
plot(zeig,’.’)
Investigate the condition number of the eigenvector matrix by computing it for ten random perturbations Q; report and comment on your results. The Matlab command eig will return eigenvalues
and the eigenvector matrix.
Notes and examples.
B = 2*rand(n) – eye(n);
[Q,R] = qr(B);
creates an n by n random matrix B with elements uniformly distributed in [−1, 1]; you can use
the QR-factorization to create a random orthogonal matrix Q. To build the matrices for your
experiments you will also need the following:
J =
0 1 0 0 0 0
0 0 1 0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 0 0 0
, S =
−2 1 0 0 0 0
1 −2 1 0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 1 −2 1 0
0 0 0 1 −2 1
0 0 0 0 1 −2
.
U = J =
0 1 0 0 0 0
0 0 1 0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 0 0 0
, LD =
−2 0 0 0 0 0
1 −2 0 0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 1 −2 0 0
0 0 0 1 −2 0
0 0 0 0 1 −2
.
J is the n by n Jordan super block, Ji,i+1 = 1 and all other elements zero. S is related to the
discrete version of the second derivative operator, S = J
T − 2I + J; U = J is the upper triangular
part of S, whereas LD = A−U = J
T −2I is the lower triangular part and diagonal of S.
mrt 1
Macm 316 Computing Assignment #4
S = -2*diag(ones(n,1)) + diag(ones(n-1,1),-1) + diag(ones(n-1,1),1);
LD = tril(S); U = triu(S,1);
Ex1 Take A = 4J. What are the exact eigenvalues and eigenvectors of A? Food for thought:
If you have the equation z
n = 0, and perturb it slightly to z
n = δ, what happens to the
solutions?
Ex2 (Limaçon) Take A = 4J + 4J
2
. What are the exact eigenvalues and eigenvectors of A?
Ex3 (Gauss-Seidel) This example is the Gauss-Seidel iteration matrix A = −(LD)
−1U (use the
Matlab command inv to compute the inverse). There is a formula for the exact eigenvalues
of A; it is known that they are real, and are all in the interval [0, 1]. Is λ = 0 an eigenvalue
of the matrix A?
Ex4 For this example we are starting out with the eigenvalues, so we definitely know what they
are – or ought to be. We take n equally spaced points on the interval [−2, 2],
tj = 4(j − 1)/(n − 1) − 2, j = 1, 2, . . . n.
We then compute the polynomial
p(z) = c0 + c1z + · · · cn−1z
n−1 + z
n =
Yn
j=1
(z − tj )
whose roots are precisely the tj . Its companion matrix is a matrix that has p(z) as its
characteristic polynomial, and therefore the tj as its eigenvalues:
A =
−cn−1 −cn−2 · · · −c1 −c0
1 0 · · · 0 0
0 1 · · · 0 0
0 0
.
.
. 0 0
0 0 · · · 1 0
.
t=4*[0:n-1]/(n-1) – 2;
p=poly(t); A=compan(p);
Ex5 Taking the same tj from Ex4, create the symmetric matrix A with tj , j = 1, 2, . . . , n as its
eigenvalues and with random (orthogonal) eigenvectors w1, . . . , wn:
B = 2*rand(n) – eye(n); [W,R]=qr(B); A = W*diag(t)*W’;
Compute the eigenvalues of A as well as those of A+δQ and A+δ(Q+QT
), for ten different
random orthogonal matrices Q. The second perturbation is symmetric, so the eigenvalues
will remain real. Do not plot your results in this case; instead, report by how much
the eigenvalues move under these perturbations. Although Ex4 and Ex5 have the same
eigenvalues, what do you believe to be the reason for obtaining very different results with
perturbations?
mrt 2
