注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

信息 灵感 创新

I? =Information,Inspiration,Innovation

 
 
 

日志

 
 
关于我

we are 5. Mathematics, Computation, Programming, Engineering, and Making fun of life.

网易考拉推荐

Jacobi方法求实对称矩阵的特征值和特征向量  

2010-12-01 16:29:01|  分类: M&M |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

%本函数采用Jacobi方法计算实对称矩阵的所有特征值和对应特征向量
%返回值D为特征值对角阵,V为对应特征向量构成的正交方阵
%即有V'*A*V=D,V'*V=I
%采用查找绝对值最大的非对角元素方法
function [D,V]=Jaco(A)
    tic;
    %检验输入是否合法
    b=size(A);
    if b(1)~=b(2)    %行列不等
        error('MATLAB:Jaco:Invalid Matrix,The Matrix input should be a Symmetry Phalanx.  See Jaco.');
    end
    n=max(b);
    for i=1:n    %非对称
        for j=1:n
            if abs(A(i,j)-A(j,i))>eps    %不能用不等号,因为机器有误差
                error('MATLAB:Jaco:Invalid Phalanx,The Phalanx input should be a Symmetric one.  See Jaco.');
            end
        end
    end
   
    %实际计算
    %初始化D,V为单位矩阵,避免多次分配空间
    %并且相乘不会造成影响
    D=eye(n);
    V=eye(n);
    %采用扫描绝对值最大A(p,q)的算法
    p=0;    %储存最大元素所在行
    q=0;    %储存最大元素所在列
    maxpq=0;    %储存绝对值最大元素
    for i=1:n-1
        for j=i+1:n
            if abs(A(i,j))>abs(maxpq)
                maxpq=A(i,j);
                p=i;
                q=j;
            end
        end
    end
    %
   
    while abs(maxpq)>eps
        maxpq=0;%务必清零,否则会死循环
        phi=atan2(2*A(p,q),A(p,p)-A(q,q))/2;
        U=eye(n);
        U(p,p)=cos(phi);
        U(q,q)=cos(phi);
        U(p,q)=-sin(phi);
        U(q,p)=sin(phi);

        D=U'*A*U;
        V=V*U;
        A=D;
       
        %改写maxpq
        for i=1:n-1
            for j=i+1:n
                if abs(A(i,j))>abs(maxpq)
                    maxpq=A(i,j);
                    p=i;
                    q=j;
                end
            end
        end
    end
    toc;
end
   

  评论这张
 
阅读(2869)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016