Posted on 2008-06-24 10:52
山泉彎延 閱讀(2480)
評論(0) 編輯 收藏 引用 所屬分類:
數(shù)值分析
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]
%通過調(diào)用malab的自帶的函數(shù)求得A的所有特征值和特征向量
%特征值保存在v中，特征向量保存的在d中,將其打印出來和我們的算法算出來的特征值進行對比
[v,d]=eig(A)
%求出行和列的大小
msize=size(A);
%取得矩陣的列數(shù)，其實行數(shù)和列數(shù)都為n
n=msize(1);
%生成n階單位陣
Q=eye(n);
%用household的方法求矩陣A的上海森伯格陣
for i=1:n-2%從第一列開始到倒數(shù)第三列
%求出每一列的最大值
d=max(abs(A(i+1:n,i)));
%規(guī)范化
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矩陣根據(jù)課本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法是否繼續(xù)進行的標志位
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的結(jié)果去替換A
A =R*Q;
%下面這個循環(huán)檢測A的精度時候足夠，去看A的次對角線各個元素的絕對值是否小于誤差限制
for w =2:n
if abs(A(w,w-1))>error
flag = 1 ;
break;%若有其中一個元素的絕對值還是大于誤差限制則還要繼續(xù)進行QR分解
end
end
%判斷的過程完畢
end
%把A打印出來
A