clc;
clear all;
%矩陣A
A=[2 ,-1,0;-1,2,-1;0,-1,2]
%取矩陣A的維數
n=max(size(A));
%叠代誤差
Eps=1E-5;
r=1;
%最大叠代次數為100
m=100;
k=1;
%小於叠代次數或叠代誤差進入計算
while r>=Eps & k<=m
p=1;
q=1;
amax=0;
for i=1:n
for j=1:n
if i~=j & abs(A(i,j))>amax
amax=abs(A(i,j));
p=i;
q=j;
end
end
end
r=amax;%計算當前叠代誤差
%以下為構造正交矩陣U
l=-A(p,q);
u=(A(p,p)-A(q,q))/2;
if u==0
w=1;
else
w=sign(u)*l/sqrt(l*l+u*u);
end
s=-w/sqrt(2*(1+sqrt(1-w*w)));
c=sqrt(1-s*s);
U=eye(n);
U(p,p)=c;
U(q,q)=c;
U(p,q)=-s;
U(q,p)=s;
%旋轉計算
A=U'*A*U%顯示每步計算A的計算結果
k=k+1;
end
if k>m
disp('A矩陣不收斂');
else
for i=1:n
D(i)=A(i,i);
end
disp('A特征值為:');
D
end
QR方法以A=[1 -1 2;-2 0 5;6 -3 6]為例不改了，自己改
構造矩陣
>>A=[1 -1 2;-2 0 5;6 -3 6]
A =
1 -1 2
-2 0 5
6 -3 6
將矩陣A變換為相似的擬上三角矩陣(即為上Hessenberg矩陣)
>>H=hess(A)
H =
1.0000 2.2136 -0.3162
6.3246 4.8000 -1.4000
0 6.6000 1.2000
對H矩陣作QR分解:
>>[Q,R]=qr(H)
Q =
-0.1562 0.2101 -0.9651
-0.9877 -0.0332 0.1526
0 0.9771 0.2127
R =
-6.4031 -5.0868 1.4322
0 6.7546 1.1526
0 0 0.3468
作50次叠代計算(具體叠代次數可依具體實驗矩陣進行)
>>for i=1:50
B=R*Q;
[Q,R]=qr(B);
end
>>R*Q
ans =
5.0000 7.4864 0.5929
-0.0000 3.0000 4.9600
0 0.0000 -1.0000
由以上結果可得到叠代計算的特征值為,可見基本QR法的叠代精度還是很高的.