【MATLAB】输入三阶矩阵A，输出矩阵W，W‘*A*W=对角矩阵

function [W,Result]=work2(A)
    %检查矩阵是否是方阵
    [rows,cols]=size(A)
    if rows ~= cols 
        disp('[work2] the matrix of input is not a square')
        return
    else if rows == cols
            if rows == 1 
                W=1
                return 
            end
            if rows ~= 3
                disp('[work2] cannot solved other matrix which dims more or less compare with three dims')
                return
            end
    end

  
    %求矩阵A的-单位特征向量构成的矩阵-和-特征值构成的对角矩阵
    [vector_matrix,diag_matrix]=eig(A)
    
    %定义U1
    U1=zeros(3,3)
    
    %U1第一列
    U1(:,1)=vector_matrix(:,1)
    
    %U2第二列
    %{
    1:  (a1,b1,c1)
    2:  (0,1,c2)
    b1*1 + c1*c2=0 解得c2=-b1/c1
    %}
    U1(:,2)=[0;1; - U1(2,1)/U1(3,1) ]
    U1(:,2)=U1(:,2)/norm(U1(:,2))
    
    %U3第三列
    U1(:,3)=null([U1(:,1)';U1(:,2)'])
    U1(:,3)=U1(:,3)/norm(U1(:,3))
    
   %
   S=U1'*A*U1
   
   %A1小矩阵
   A1=S(2:3,2:3)
    
   %求矩阵A1的-单位特征向量构成的矩阵-和-特征值构成的对角矩阵
   [vector_matrix1,diag_matrix1]=eig(A1)
   
   %定义V
   V=zeros(2,2)
   
   %V第一列
   V(:,1)=vector_matrix1(:,1)
   
   %V第二列
   V(:,2)=null(V(:,1)')  %单个方程求解
   V(:,2)=V(:,2)/norm(V(:,2))
   
   %第一行乘上V
   tmp=S(1,2:3)*V
   
   %
   S1=V'*A1*V
   
   %定义U2
   U2=zeros(3,3)
   U2(1,1)=1
   U2(2:3,2:3)=V
   
   %
   W=U1*U2
   
   %定义Result
   Result=zeros(3,3)
   Result(1,1)=S(1,1)
   Result(1,2:3)=tmp
   Result(2:3,2:3)=S1
end