% Script 3 for Week 4
% The Power Method
% We want to find evecs and evals for matrix A.
% These are numbers L and nonzero vectors x for which A*x=L*x.
% So A has the simplest possible effect on evecs: it acts by numerical multiplication by a scalar.
% L is an eval when (L*eye(n)-A)*x=0 has a nontrivial solution which happens exactly when det(L*eye(n)-A)=0.
A=[40 -20 0
-20 40 -20
0 -20 40]
f=@(L) det(L*eye(3)-A)
fplot(f,[-100,100])
fplot(f,[0,90])
hold on
plot([0 100],[0,0])
hold off
% By inspection you can see that A has three eigenvalues located near 12, 39 and 68.
% eigpics will pick a random seed vector and successively magnify (here 6 times) the fraction of x that is a multiple of the evec for the largest eval.
x=eigpics(A,6)
x=x/x(1)
y=A*x
y=y/y(1)
norm(x-y)
% x is pretty close to an evec--nearly 99% evec for an eval near 68. The rapidity of convergence depends on the hugeness of the biggest eval in comparison to the second hugest and the fraction of the initial random seed in the eigenvector-direction for that largest eval.
% The graph shows that there is another eval near 11.5. We would like an evec for that one too.
% Note that for any real c the matrices A-c*eye(3) and inv(A-c*eye(3)) have the same evecs as does A.
% If x,L is an evec-eval pair for A then x,L-c will be for A-c*eye(3)
and x,1/(L-c) will be for inv(A-c*eye(3)).
% If c is close to L then eval 1/(L-c) for inv(A-c*eye(3)) will be gigantic so eradication of OTHER evec components will be very rapid.
% Going back to the same matrix A from above
B=inv(A-11.5*eye(3));
% Run eigpics on B.
x=eigpics(B)
x=x/x(1)
y=A*x
eigenval=y(1)
y=y/y(1)
norm(x-y)
% The easy way to find evals/vecs???
e=eig(A)
[E,e]=eig(A)
% Another example:
P=[1 1 3; 4 5 1; 6 1 0]
D=diag([-2 2 -3])
H=P*D*inv(P)
% (I am using P and D to produce the matrix H with known evals and evecs.)
x=eigpics(H)
[E,e]=eig(H)
% Complex e-vals of real matrices COME IN CONJUGATE PAIRS. Corresponding e-vecs also have complex entries, and MAY BE CHOSEN TO BE CONJUGATES OF EACH OTHER.
% In many of the applications of this stuff answers must be real, so coordinates conspire to cancel out the complex parts of vectors and measurable values.
P=[1 1 3; 4 5 1; 6 1 0]
W=P*[-2 0 0;
0 2 -3;
0 3 2]*inv(P)
% (I am using P to produce the matrix W with known complex evals and evecs.)
[E,e]=eig(W)
% Here, the largest magnitude eval is complex so the direct power method doesn't work:
[x,rotation,rotangs]=eigpics(W)
2*pi./rotangs
[x,rotation,rotangs]=eigpics(W,1000);
figure(7)
plot(rotangs)
% As you can see the x vector never seems to "settle down", though the vectors become (nearly) co-planar very quickly. The angle between succesive vectors seems to be an irrational multiple of pi.
% But .... WHY doesn't it settle down? A generic real seed vector will be a linear combination of the three eigenvectors. Since the real eigenvalue is smallest that component gets crushed flat by the iteration, and what remains is a combination of the two complex evecs. But since the combination must be real, the coefficients of that combination must be conjugates. Anyway, we have
x=a*v+(conjugate a)*(conjugate v).
% (Note: it is easy to show that if v=v_r+i*v_c breaks vector v into its real and imaginary parts and if a=s+i*t then x=2*s*v_r-2*t*v_c.)
% So A*x=(lambda)*a*v+(conjugate lambda)*(conjugate a)*(conjugate v) and generally
% A^n*x=(lambda)^n*a*v+(conjugate lambda)^n*(conjugate a)*(conjugate v)
% We are bouncing around in the plane of v and conjugate v. But there is a pattern.
% Lambda rotates each ENTRY of v by a certain angle in the complex plane and (conjugate lambda) rotates each entry of (conjugate v) by the negative of that angle. If the angle of lambda above is NOT a rational multiple of pi the rotations will never return any of the evec entries to their original angles. Their combination, which will be a real vector, changes every time it is hit by A and never repeats.
% But IS the angle ACTUALLY an irrational multiple of pi? It is fairly easy to show (using De Moivre's Theorem) that if r=p/q is rational then cos(r*pi) is an algebraic number: that is, the root of a polynomial equation of degree 2*q.
% Specifically, let c=cos(r*pi) and s=sin(r*pi).
% exp(p*pi/q)^(2*p)=(c+i*s)^(2*p)=cos(2*p*pi)+i*sin(2*p*pi)=1.
% So the complex part of the expansion of (c+i*s)^(2*p) is zero and the real part is
% 1 = sum of nchoosek(2*q,k)*c^k(-1)^((2*q-k)/2)*(1-c^2)^((2*q-k)/2)
% where the sum is over EVEN k only and nchoosek(2*q,k) is the binomial coefficient.
% But WHICH algebraic numbers c correspond to rational multiples of pi?
% A partial answer is given in 2006 in Juan L. Varona's paper "Rational Values of the Arccosine Function." There, Verona proves that if x is rational with 0