带wilkinson位移的隐式对称qr算法matlab代码

以下是带 Wilkinson 位移的隐式对称 QR 算法的 MATLAB 代码：

function [Q, H] = symmqr(A, tol)
% SYMMQR  Implicit symmetric QR algorithm with Wilkinson shift for eigenvalue computation.
%    [Q, H] = SYMMQR(A, TOL) computes the eigenvalues of the input matrix A
%    using the implicit symmetric QR algorithm with Wilkinson shift,
%    where TOL is the tolerance for convergence.
%    Q is the orthogonal matrix satisfying A = Q*H*Q', where H is the
%    tridiagonal matrix with the eigenvalues on the diagonal.

n = size(A, 1);
H = hessenberg(A);
Q = eye(n);

while norm(H(n, 1:n-1)) > tol
    % Compute Wilkinson shift
    mu = H(n, n);
    a = H(n-1, n-1);
    b = H(n-1, n);
    c = H(n, n-1);
    delta = (a - mu)*(a - mu) + b*b;

    if abs(delta) < eps
        t = 1;
    else
        t = (a - mu + sign(a - mu)*sqrt(delta)) / b;
    end

    % QR step with Wilkinson shift
    for k = 1:n-1
        if k == n-1 || abs(H(k+1, k)) <= tol
            % Single QR iteration
            x = [H(k, k) - t; H(k+1, k)];
            if abs(x(1)) < abs(x(2))
                c = x(1); s = x(2);
            else
                tau = x(2) / x(1);
                c = 1 / sqrt(1 + tau*tau);
                s = tau*c;
            end

            % Apply Givens rotation to H and Q
            G = [c, s; -s, c];
            H(k:k+1, k:n) = G'*H(k:k+1, k:n);
            H(1:k+1, k:k+1) = H(1:k+1, k:k+1)*G;
            Q(:, k:k+1) = Q(:, k:k+1)*G;
        else
            % Double QR iteration
            x = [H(k, k) - t; H(k+1, k); H(k+2, k)];
            if abs(x(2)) < abs(x(3))
                tau = x(2) / x(3);
                s = sqrt(1 / (1 + tau*tau));
                c = tau*s;
            else
                tau = x(3) / x(2);
                c = sqrt(1 / (1 + tau*tau));
                s = tau*c;
            end

            % Apply Givens rotation to eliminate subdiagonal element
            G = [c, s, 0; -s, c, 0; 0, 0, 1];
            H(k:k+2, k:n) = G'*H(k:k+2, k:n);
            H(1:k+2, k:k+2) = H(1:k+2, k:k+2)*G;
            Q(:, k:k+2) = Q(:, k:k+2)*G;
        end
    end
end

end

function H = hessenberg(A)
% Hessenberg reduction of matrix A using Householder transformations.
n = size(A, 1);

for k = 1:n-2
    % Compute the Householder reflection vector
    x = A(k+1:n, k);
    e1 = zeros(length(x), 1);
    e1(1) = 1;
    v = sign(x(1))*norm(x)*e1 + x;
    v = v / norm(v);

    % Apply the Householder reflection to A
    A(k+1:n, k:n) = A(k+1:n, k:n) - 2*v*(v'*A(k+1:n, k:n));
    A(1:n, k+1:n) = A(1:n, k+1:n) - 2*(A(1:n, k+1:n)*v)*v';
end

H = triu(A, -1);
end

其中，`A` 是输入的矩阵，`tol` 是算法收敛的阈值，`Q` 是正交矩阵，满足 $A=QHQ^T$，`H` 是带有特征值的对称三对角矩阵。该算法使用 Hessenberg 约化和隐式对称 QR 算法，以及 Wilkinson 位移来计算特征值。具体实现细节可以参考代码中的注释。