#% FILE anton6.mpl
interface(echo=2);

# AR6 Section 6.1 Example 2 p269
with(linalg,matrix,eigenvals):
A2:= matrix([ [ 3, 2],
                [-1, 0] ]):
  A:=A2:
eigenvals(A2);

# AR6 Section 6.1 Example 3 p269; AR7 Section 7.1 p358
A3:= matrix([ [-2,-1],
                [ 5, 2] ]):
eigenvals(A3);
B:= matrix([ [0,-1],
               [1, 0] ]):
eigenvals(B);


# {EXERCISE.} Find the eigenvalues of the following 
#  three symmetric tridiagonal matrices:
  C1:= matrix([[1,1,0],[1,0,1],[0,1,1]]);
  C2:= matrix([[0,1,0],[1,0,1],[0,1,0]]);          
  C3:= matrix([[1,1,0],[1,0,1],[0,1,0]]);

evalf(Eigenvals(C3));

# Continuing the preceding example:
with(linalg,diag,nullspace):
A:=A2: Id:= diag(1,1):                W1:= nullspace( evalm(A - Id) );
W_2:= nullspace( evalm(A - 2*Id) );

# For example,  W1  is a basis for  nullSpace(A-Id),
# and hence a set of eigenvectors.
# To CHECK that this set of vectors satisfies the equation  A x = x,  find
evalm(A &* W1[1]);   # The product of A with the first (and only) vector in W1.

# You should CHECK the other eigenvectors for yourself.

# Example 5 - AR6 Section 6.1 p271; Section 7.1 p359
with(linalg,matrix,eigenvals,eigenvects,augment,equal,inverse,diag):
A:= matrix([ [0,0,-2],
               [1,2, 1],
               [1,0, 3] ]):
evects:= eigenvects(A);

# The eigenvectors are the elements of the sets  evects[1][3] and evects[2][3]:
if (nops(evects[1][3])=1) then
  evals:= evects[1][1], evects[2][1], evects[2][1] :
  P:=augment( evects[1][3][1], evects[2][3][1], evects[2][3][2] );
else
  evals:= evects[1][1], evects[1][1], evects[2][1] :
  P:=augment( evects[1][3][1], evects[1][3][2], evects[2][3][1] );
fi;

equal( evalm( inverse(P)&*A&*P ), diag(evals) );

# Not all matrices are diagonalizable.  For Example 
B:= matrix([ [a,1], [0,a] ]):
evals:= eigenvals(B);        # eigenvects will not work on this matrix

nullspace( B-a*diag(1,1) );
# Since the dimension of this eigenspace is only 1,   B is not diagonalizable.


with(linalg):
u:= vector([I,1+I]):            v:= vector([-2*I,1+I]):
simplify(dotprod(u,v));
innerprod(u,v);


# AR6 Section 6.3 Example 2 p288; AR7 Section 7.3 Example 1 p377
with(linalg,matrix,eigenvals,eigenvects,dotprod,
              GramSchmidt,norm,augment,orthog,diag,equal):
# Find an orthogonal matrix Q which diagonalizes the matrix A where
A:= matrix([ [4, 2, 2],
               [2, 4, 2],
               [2, 2, 4] ]):
# Use the same ideas as in an earlier example:
evects:= eigenvects(A);
if (nops(evects[1][3])=1) then 
  evals:= evects[1][1], evects[2][1], evects[2][1] :
  u1:= evects[1][3][1]: u2:= evects[2][3][1]: u3:= evects[2][3][2]:
else
  evals:= evects[2][1], evects[1][1], evects[1][1] :
  u1:= evects[2][3][1]: u2:= evects[1][3][1]: u3:= evects[1][3][2]:
fi;
print( dotprod(u1,u2), dotprod(u1,u3), dotprod(u2,u3) );

# Since  u2, u3  are not orthogonal, use Gram-Schmidt:
v1:= evalm( u1/norm(u1,2) ):          GS:= GramSchmidt( [u2,u3] ):
v2:= evalm( GS[1]/norm(GS[1],2) ):    v3:= evalm( GS[2]/norm(GS[2],2) ):
Q:= augment( v1,v2,v3 );
orthog(Q);
equal( evalm( transpose(Q)&*A&*Q ), diag(evals) );


with(linalg):
# Show that Q reflects vectors orthogonal to some plane where
Q:= matrix([ [4/9,-7/9,-4/9],
             [1/9,-4/9, 8/9],
             [8/9, 4/9, 1/9] ]):
print( orthog(Q), det(Q) );

# Thus  Q reflects vectors orthogonal to a plane  P  where
u:= nullspace( Q+diag(1,1,1) )[1]:   # u is an eigenvector associated with -1
P:= dotprod( u,vector([x,y,z]) ) = 0;

# It also rotates the plane  P  through an angle  theta  where
cos_theta:= (trace(Q)+1)/2;


with(linalg,band,eigenvals):
for  n from 2 to 8 do  A.n:=band([-1,2,-1],n): eigenvals(A.n); od;



# TIDY THE Maple MESS BELOW

with(linalg,eigenvals,eigenvects,band,charmat,kernel,scalarmul,norm):
n:=3:
A.n:=band([-1,2,-1],n);
# Eigenvals and their eigenvectors can be found numerically
readlib(forget): forget(Eigenvals):
evalf(Eigenvals(A3,vec3));
print(vec3);
# eigenvals and their eigenvects can sometimes be found exactly 
evects:=eigenvects(A.n);
# Maple has no real problem with the algebraic numbers here
x:='x': solve(x^2 - 4*x +2 =0, x); # for n=3
tidy:= g -> op(map(op,g));
evecta:=[tidy([allvalues(evects[1][3])]),tidy([allvalues(evects[2][3])])];
evectf:=map(evalf,evecta);
nmze:= vec -> scalarmul(vec,1/norm(vec,2));
evectnf:=map(nmze,evectf);
# compare the two approaches: they agree
print(vec3);
# Another way to get eigenvects exactly without so many allvalues is
charmat(A3,2+sqrt(2)); # for n=3
kernel(");

# ADD THE CAUTION ABOUT readlib(forget) STUFF

# Synge and Griffith S17.3 p500
########################################################################
# n-masses with equal springs joining them.
#
# p501 xkj is the k-th element of j-th eigenvector
n:='n': xkj:= sin(k*j*Pi/(n+1)); # really *sqrt(2/(n+1)) to normalise
# let A=band([-1,2,-1],n) and show xj is an eigenvector:
# find corresponding eigenvalue.
Ax:=expand(-subs(k=i-1,xkj) +2*subs(k=i,xkj) - subs(k=i+1,xkj));
#                        i j Pi         j Pi            i j Pi
#               - 2 sin(--------) cos(-------) + 2 sin(--------)
#                         n + 1        n + 1             n + 1
#
evalj:=factor(Ax)/subs(k=i,xkj);
#                                       j Pi
#                             - 2 cos(-------) + 2
#                                      n + 1
#
# To fix n and plot out the eigenvectors use
# nPts:=7; for jp from 1 to np do  # if all 8 eigenvects wanted
nPts:= 7: jp:=4:
ptlist:=[0,0]:
for kp from 1 to np do
  ptlist:=[op(ptlist),kp,subs(n=nPts,k=kp,j=jp,xkj)]:
  od:
ptlist:=[op(ptlist),nPts+1,0]:
interface(plotdevice=postscript,plotoutput=EVplotfile);
plot(ptlist,style=LINE);

interface(echo=1);
