Posted on 2008-06-24 10:52
山泉弯延 阅读(2448)
评论(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 ;
R =A;
Q = 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
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