eig_jacobi.m

求对称矩阵全部特征值的Jacobi方法

function [D,R]=eig_Jacobi(A, ep)
% 求对称矩阵全部特征值的Jacobi方法，其中
%    A      --- 矩阵
%    ep     --- 精度要求
%    D      --- 对角线上的值为特征值
%    R      --- 对应特征值的特征向量

if nargin <2  ep=1e-5;    end
n=length(A); R=eye(n);
while 1
    % 求最大元素
    Amax=0;
    for l=1:n-1
        for k=l+1:n
            if abs(A(l,k))>Amax
                Amax=abs(A(l,k)); i=l; j=k;
            end
        end
    end
    % 终止准则
    if Amax<ep  break;  end
    % 计算三角函数
    d=(A(i,i)-A(j,j))/(2*A(i,j));
    if abs(d)<1e-10
        t=1;
    else
        t=sign(d)/(abs(d)+sqrt(d^2+1)); 
    end
    c=1/sqrt(t^2+1); s=c*t;
    % 旋转计算
    for l=1:n
        if l==i 
            Aii=A(i,i)*c^2+A(j,j)*s^2+2*A(i,j)*s*c;
            Ajj=A(i,i)*s^2+A(j,j)*c^2-2*A(i,j)*s*c;
            A(i,j)=(A(j,j)-A(i,i))*s*c+A(i,j)*(c^2-s^2);
            A(j,i)=A(i,j); A(i,i)=Aii; A(j,j)=Ajj;
        elseif l~=j
            Ail=A(i,l)*c+A(j,l)*s;
            Ajl=-A(i,l)*s+A(j,l)*c;
            A(i,l)=Ail; A(l,i)=Ail;
            A(j,l)=Ajl; A(l,j)=Ajl;
        end
        Rli=R(l,i)*c+R(l,j)*s;
        Rlj=-R(l,i)*s+R(l,j)*c;
        R(l,i)=Rli; R(l,j)=Rlj;
    end
end
D=diag(diag(A));