newplan

阿基米德在洗澡時發現浮力原理,高興得來不及穿㆖褲子,跑到街㆖大喊:Eureka(我找到了)。
posts - 39, comments - 26, trackbacks - 0, articles - 4
  C++博客 :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理

数值分析方阵的QR分解

Posted on 2008-06-24 10:52 山泉弯延 阅读(2051) 评论(0)  编辑 收藏 引用 所属分类: 数值分析
function[]=iqr()
% 实验名称:方阵的QR分解
% 实验描述:先将方阵化为上海申博格阵,再用QR分解法求上海申博格阵的特征值,则所得到的特征值也是方阵的特征值
% 作者:newplan
% 实验完成日期:6月10号
%下面的A为测试三阶的方阵
A
=[5,-3,2;6,-4,4;4,-4,5]
%下面的A为测试四阶的方阵
%A 
= [1 2 1 2;2 2 -1 1;1 -1 1 1;2 1 1 1]
%通过调用malab的自带的函数求得A的所有特征值和特征向量
%特征值保存在v中,特征向量保存的在d中,将其打印出来和我们的算法算出来的特征值进行对比
[v,d]
=eig(A)
%求出行和列的大小
msize
=size(A);
%取得矩阵的列数,其实行数和列数都为n
n
=msize(1);
%生成n阶单位阵
Q
=eye(n);
%用household的方法求矩阵A的上海森伯格阵
for i=1:n-2%从第一列开始到倒数第三列 
    %求出每一列的最大值
    d
=max(abs(A(i+1:n,i)));
    %规范化
    U(i
+1:n,i)=A(i+1:n,i)/d;
    delta
=U(i+1,i)*norm(U(i+1:n,i))/abs(U(i+1,i));
    U(i
+1,i)=U(i+1,i)+delta;
    beta 
= delta*U(i+1,i);
    %求出R矩阵根据课本316P例题三 
    R 
= eye(n-i,n-i)-inv(beta)*U(i+1:n,i)*U(i+1:n,i)';
    u=eye(n,n);
    
for j =i+1:n
        
for k =i+1:n
            u(j,k)
=R(j-i,k-i);
        
end
    
end
    A
=u*A*u;%生成新的A=u×A×u
end
%error为我们设定的误差限制
error = 0.0000001;
%flag为判断QR法是否继续进行的标志位
flag 
=1;
while flag==1
flag 
=0 ;
=A;
= eye(n,n);
%按照QR分解法求出cos,
sin 然后计算V,最终得到R和Q
for i=1:n-1
  r 
= norm(R(i:i+1,i));
  icos
=R(i,i)/r;
  isin
=R(i+1,i)/r;
  v
=eye(n,n);
  v(i,i)
=icos;
  v(i
+1,i+1)=icos;
  v(i,i
+1)=isin;
  v(i
+1,i)=-isin;
  R
=v*R;
  Q
=Q*v';
end
%用R
*Q的结果去替换A
=R*Q;
%下面这个循环检测A的精度时候足够,去看A的次对角线各个元素的绝对值是否小于误差限制
for w =2:n
     
if abs(A(w,w-1))>error
     flag 
= 1 ;
     break;%若有其中一个元素的绝对值还是大于误差限制则还要继续进行QR分解
    
end   
end
%判断的过程完毕
end
%把A打印出来
A




只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理