QLmethod.m

% QL方法
% 注释的地方需要修改

function QLmethod
A=[6  2  3  4  15;
   2  8  6  8  24;
   3  6  9  9  27;
   4  8  9 16  37;
  15 24 27 37 103];         % 所求矩阵
A0=householder(A);
F=1;
while F>1e-12
    [Q,L]=QL(A0);
    A0=L*Q;
    F=norm(A0-diag(diag(A0)));
end
D=diag(A0);
[V]=modes(A,D);
output(A,V,D)

function A=householder(A)
n=length(A);
for k=1:n-2
    u=A(k+1:end,k);
    u(1)=u(1)+sign(u(1))*norm(u,2);
    beta=norm(u,2)^2/2;
    U=eye(n);
    U(k+1:end,k+1:end)=eye(length(u))-u*u'/beta;
    A=U*A*U;
end   

function [Q,L]=QL(A)
L=A;
n=length(A);
Q=eye(n);
for k=n:-1:2
    r=sqrt(L(k,k)^2+L(k-1,k)^2);
    c=L(k,k)/r;
    s=-L(k-1,k)/r;
    R=eye(n);
    R(k-1,k-1)=c; R(k-1,k)=-s;
    R(k,k-1)=s;   R(k,k)=c;
    L=R'*L;
    Q=Q*R;
end

function [V]=modes(A,D)
V=[];
for k=1:length(D)
    B=A-D(k)*eye(size(A));
    B22=B(2:end,2:end);
    B21=B(2:end,1);
    x=[1;-B22\B21];
    V=[V,x/norm(x)];
end

function rho=shift(A)
d1=A(1,1); d2=A(2,2); e1=A(1,2);
delta=(d2-d1)/2
rho=d1+delta-sign(delta)*sqrt(delta^2+e1^2)

function output(A,V,D)
n=length(A);
for i = 1:n
    fprintf('lambda%1d = %.8f\nv%1d = ',i,D(i),i); % 结果有效数字位数为.和f之间的数
    for j = 1:n
        fprintf('%10.8f ',V(j,i));                   % 结果有效数字位数为.和f之间的数
    end
    fprintf('\n\n');
end