【matlab】 QR分解 求矩陣的特征值

"QR_H.m"

function [Q,R] = QR_tao(A)
%輸入矩陣A
%輸出正交矩陣Q和上三角矩陣R
[n,n]=size(A);
E = eye(n);
X = zeros(n,1);
R = zeros(n);

P1 = E;
for k=1:n-1
    s = -sign(A(k,k))*norm(A(k:n,k));
    R(k,k) = -s;
    if k == 1
        w = [A(1,1)+s,A(2:n,k)']';
    else
        w = [zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';
        R(1:k-1,k) = A(1:k-1,k);
    end
    if norm(w)~=0
        w = w/norm(w);
    end
    P = E-2*w*w';
    A = P*A;
    P1 = P*P1;
    R(1:n,n) = A(1:n,n);
end

 

之后根據算法:

An = Q1*R1;

An+1 =  R1*Q1

重復迭代即可。

"QR.m"

%輸入 矩陣A 和迭代次數 it_max
%輸出 最后對角線上元素為特征值的矩陣
function [Q] = QR(A,it_max)
A1 = A;
for N=1:it_max
[Q1,R1] = QR_tao(A1);
A2 = R1*Q1;
A1 = A2;
end
Q=A1

測試： 計算一個矩陣的特征值:

A = [13,-3,-1,0;
    -3,6,0,-2;
    -1,0,10,-1;
    0,-2,-1,7;
];
[Q] = QR(A,50)
eig(A)

 

最后結果: