disp('PD decomposition of a matrix A')
a=[3 -2 4;-2 6 2;4 2 3]
disp(a,'A=')
disp('Eigenvalues of A are')
eig=spec(a)
disp(eig)
disp(eig(2,1),'for lambda =')
disp('A-(lambda)I=')
b=a-eig(2,1)*eye(3,3)
disp(b)
disp('To find eigenvector, form an augmented matrix')
c=[b [0;0;0]]
disp(c)
disp('performing row operations')
c(2,:)=c(2,:)-(c(2,1)/c(1,1))*c(1,:)
c(3,:)=c(3,:)-(c(3,1)/c(1,1))*c(1,:)
disp(c)
disp('With x2 and x3 as free variables, we get two vectors.')
disp('x1=-.5x2+x3')
disp('Thus, the two vectors are')
v1=[-1;2;0]
v2=[1;0;1]
disp(v2,v1)
disp('Orthogonalizing v1 and v2')
disp('Let x1=v1')
disp('x2=v2-((v2.v1)/(v1.v1))*v1')
x1=v1
a1=v2'*v1
a2=v1'*v1
x2=v2-(a1/a2)*v1
x1=x1/(sqrt(x1(1,1)^2+x1(2,1)^2+x1(3,1)^2))
x1=x2/(sqrt(x2(1,1)^2+x2(2,1)^2+x2(3,1)^2))
disp('An orthonormal basis is:')
disp(x2,x1)
disp(eig(1,1),'for lambda=')
disp('A-(lambda)I=')
b=a-eig(1,1)*eye(3,3)
disp(b)
disp('To find eigenvector, form an augmented matrix')
c=[b [0;0;0]]
disp(c)
disp('performing row operations')
c(2,:)=c(2,:)-(c(2,1)/c(1,1))*c(1,:)
c(3,:)=c(3,:)-(c(3,1)/c(1,1))*c(1,:)
disp(c)
c(3,:)=c(3,:)-(c(3,2)/c(2,2))*c(2,:)
disp(c)
c(1,:)=c(1,:)/c(1,1)
c(2,:)=c(2,:)/c(2,2)
disp(c)
c(1,:)=c(1,:)-(c(1,2)/c(2,2))*c(2,:)
disp(c)
disp('With x3 as free variable')
disp('x1=x3 and x2=-.5x3')
disp('Thus a basis for the eigenspace is:')
v3=[1;-.5;1]
disp(v3)
disp('upon normalizing')
v3=v3/(sqrt(v3(1,1)^2+v3(2,1)^2+v3(3,1)^2))
disp(v3)
disp('Thus, matrix P=')
disp([x1 x2 v3])
disp('Corresponding matrix D=')
disp([eig(2,1) 0 0;0 eig(3,1) 0;0 0 eig(1,1)])