eig2.m

实现共轭对称矩阵的特征分解算法

clear
clc
p=4; %共轭对称方阵维数
R = [2.562847+j*0,1.880512-j*0.007008409,2.095236-j*0.01132958,2.065702-j*0.1087133;1.880512+j*0.007008181,2.563189-j*0,2.076373+j*0.05680992,1.925068-j*0.04684043;2.095236+j*0.01132934,2.076373-j*0.05681005,2.919492-j*0,2.320824-j*0.07811832;2.065702+j*0.108713,1.925068+j*0.04684028,2.320824+j*0.07811818,2.818176-j*0]; %待分解矩阵

Rx = R;
U0 = eye(p);
for k=1:(p-2)
    cc = Rx(k+1:p,k); %/max(abs(Rx(k+1:p,k)));
    sigma = sqrt(sum(abs(cc).^2));
    u = cc;
    u(1) = u(1) + sigma;
    b = sigma*(cc(1)+sigma);
    r = eye(p-k)-u*u'/b;
    U = [eye(k) zeros(k,size(r,1));zeros(size(r,1),k) r];
    Rx = U*Rx*U';
    U0 = U*U0;
end

eps=0.0001;
for t=1:p-1
    upper_sum = (sum(sum(abs(Rx).^2))-sum(abs(diag(Rx)).^2))/2;
    if upper_sum>eps
       ss=Rx(p-t+1,p-t+1);
       Rx(1:p-t+1,1:p-t+1)=Rx(1:p-t+1,1:p-t+1)-ss*eye(p-t+1);
       for k=1:p-t
           if abs(Rx(k+1,k))>eps
              g=sqrt(abs(Rx(k,k))^2+abs(Rx(k+1,k))^2); %g=sigma
              c=Rx(k,k)/g;
              s=Rx(k+1,k)/g;
              pp=[c,s';-s,c];
              P=[eye(k-1),zeros(k-1,2),zeros(k-1,p-k-1);zeros(2,k-1),pp,zeros(2,p-k-1);zeros(p-k-1,k-1),zeros(p-k-1,2),eye(p-k-1)];
              U0=P*U0;
              Rx=P*Rx*P';
           end
       end
       Rx(1:p-t+1,1:p-t+1)=Rx(1:p-t+1,1:p-t+1)+ss*eye(p-t+1);
    end
end
disp('特征向量为')
U0
disp('对应特征值')
diag(Rx)