jac.m

雅克比方法求矩阵的特征值和特征向量

function[aa,v]=jac(a)
%此程序用jacobi方法求实对称矩阵的全部特征值和特征向量
%输入x：n x n矩阵
%----------------------------------------------------%
    n=length(a);
    aa=zeros(1,n);
    u=zeros(n,n);
    
    l=0;
    %v是n阶单位矩阵
    v=eye(n,n);
    while (1)
        l=l+1;
        %求绝对值最大的矩阵元素fm,第p行，第q列
        fm=0.0;
        for i=1:n
            for j=1:n
                if i==j 
                    continue;
                else
                d=abs(a(i,j));
                if d>fm
                    fm=d;
                    p=i;
                    q=j;
                end
                end
            end
        end
        if fm<0.0001
            break;
        end
        x=2*a(p,q)*(sign(a(p,p)-a(q,q)));
        y=abs(a(p,p)-a(q,q));
        if y~=0
           cn=sqrt((1+y/sqrt(x^2+y^2))/2.0);
        else
            cn=1/sqrt(2);
        end
        if y~=0
            sn=x/sqrt(x^2+y^2)/(2*cn);
        elseif y==0 & a(p,q)>0
            sn=1/sqrt(2);
        elseif y==0 & a(p,q)<0
            sn=-1/sqrt(2);
        end
        %更新a
        a_t=a;
        for i=1:n
            if i~=p & i~=q
                a(i,p)=a_t(i,p)*cn+a_t(i,q)*sn;
                a(i,q)=-a_t(i,p)*sn+a_t(i,q)*cn;
                a(p,i)=a(i,p);
                a(q,i)=a(i,q);
            end
        end

        a(p,p)=a_t(p,p)*cn*cn+2*a_t(p,q)*cn*sn+a_t(q,q)*sn*sn;
        a(q,q)=a_t(p,p)*sn*sn-2*a_t(p,q)*sn*cn+a_t(q,q)*cn*cn;
        a(p,q)=0.0;
        a(q,p)=a(p,q);
        %更新v
        v_t=v;
        for i=1:n
            v(i,p)=v_t(i,p)*cn+v_t(i,q)*sn;
            v(i,q)=-v_t(i,p)*sn+v_t(i,q)*cn;
        end
    end
    for i=1:n
        aa(1,i)=a(i,i);
    end
end