jocobirot.m

本程序是根据jacobi旋转法求实对称矩阵的全部特征值和特征向量

function [EigValMat,EigVecMat]=JocobiRot(A,eps)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 本程序是根据jacobi旋转法求实对称矩阵的全部特征值和特征向量
%% 输入变量:A为对称实矩阵,eps为允许误差.
%% 输出变量:EigValMat为A的特征值,EigVecMat为特征向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(A);n=n(1);%% n为矩阵A的阶数
P=eye(n);%% P为旋转矩阵,赋初值为单位阵
trc=1;%% trc为矩阵A的非对角元素的平方和,赋初值为1;
num=0;%% num设置为累加器,记录迭代的次数
while abs(trc)>=eps%% 进行正交变换A=PAP'将A的两个绝对值最大的非对角元素化零,直到所有非对角元素的平方和小于给定的eps,则结束循环.
    MaxMes=FinMax(A);%% 寻找绝对值最大的非对角元素的位置及所有非对角元素的平方和
    l=MaxMes(1);%% l为绝对值最大的非对角元素的行号
    s=MaxMes(2);%% s为绝对值最大的非对角元素的列号
    trc=MaxMes(3);%% trc为矩阵A的非对角元素的平方和
    
    RotMat=ComRotMat(A,l,s);%% 计算此次旋转的平面旋转矩阵RotMat
    A=RotMat*A*RotMat';%% 对当前A进行一次旋转变换将A的两个绝对值最大的非对角元素化零,并仍记为A
    P=RotMat*P;%% 记录到目前为止所有旋转矩阵的乘积
    num=num+1;%% 记录已经进行旋转的次数
end
EigValMat=A;
EigVecMat=P;
num  

function MaxMes=FinMax(A) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%对上三角进行扫描，寻找矩阵A的绝对值最大的非对角元素的位置及所有非对角元素的平方和
%%输入量:实对称矩阵A
%%输出量:MaxMes记录矩阵A的绝对值最大的非对角元素的行号列号及所有非对角元素的平方和
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(A);n=n(1);
trc=0;
l=1;s=2;
max=abs(A(l,s));
for i=1:n-1
    for j=i+1:n
        if max<abs(A(i,j))
            max=abs(A(i,j));
            l=i;
            s=j;
        end
        trc=trc+2*A(i,j)^2;
    end
end
MaxMes=[l,s,trc]';


function RotMat=ComRotMat(A,l,s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%计算平面旋转矩阵
%%输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s
%%输出量:旋转的平面旋转矩阵RotMat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(A);
n=n(1);
RotMat=eye(n);
if A(l,l)==A(s,s)
    if A(l,s)>0
        cos1=1/sqrt(2);
        sin1=-1/sqrt(2);
    else 
        cos1=1/sqrt(2);
        sin1=1/sqrt(2);
    end
else 
    t=-2*A(l,s)/(A(l,l)-A(s,s));
    cos1=sqrt((1/sqrt(1+t^2)+1)/2);
    sin1=t/(2*sqrt(1+t^2)*cos1);
end

RotMat(l,l)=cos1;
RotMat(l,s)=sin1;
RotMat(s,l)=-sin1;
RotMat(s,s)=cos1;