带Wilkinson位移的隐式对称QR算法matlab代码

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

function [eigenvalues, eigenvectors] = wilkinson_shift_hessenberg(A, tol)
% Compute eigenvalues and eigenvectors of a symmetric matrix A using the
% implicit symmetric QR algorithm with Wilkinson shift.
% Input: A - symmetric matrix
%        tol - tolerance for convergence
% Output: eigenvalues - vector of eigenvalues
%         eigenvectors - matrix of eigenvectors

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

while true
    % Compute the Wilkinson shift
    mu = H(n, n);
    a = H(n-1, n-1);
    b = H(n, n-1);
    c = H(n-1, n);
    delta = (a - mu)^2 + b^2 - c*(a + mu);
    if delta < 0
        sgn = -1;
    else
        sgn = 1;
    end
    if abs(b) > abs(a-mu)
        t = -c / (b + sgn*sqrt(delta));
    else
        t = (a-mu) / (b + sgn*sqrt(delta));
    end
    x = H(1,1) - t;
    y = H(2,1);
    
    % QR iteration with Wilkinson shift
    for k = 1:n-1
        [Q, R] = qr(H(1:n-k+1, 1:n-k+1) - x*eye(n-k+1));
        H(1:n-k+1, 1:n-k+1) = R*Q + x*eye(n-k+1);
        H(1:n-k+2, n-k+2) = 0;
        if k > 1
            H(k, k-1) = 0;
        end
        eigenvectors(:, 1:n-k+1) = eigenvectors(:, 1:n-k+1)*Q;
        x = H(k+1, k);
        y = H(k+2, k);
        if abs(y) < tol*(abs(x) + abs(H(k+2, k+1)))
            H(k+2, k) = 0;
            break;
        end
        % Compute the Givens rotation to zero out y
        [c, s] = givens_rotation(x, y);
        G = eye(n);
        G(k:k+1, k:k+1) = [c, s; -s, c];
        H(1:n-k+2, k:k+1) = H(1:n-k+2, k:k+1)*G;
        H(k:k+1, 1:n-k+2) = G'*H(k:k+1, 1:n-k+2);
        eigenvectors(:, k:k+1) = eigenvectors(:, k:k+1)*G;
    end
    
    % Check for convergence
    if k == n-1
        eigenvalues = diag(H);
        break;
    end
end
end

function [c, s] = givens_rotation(x, y)
% Compute the Givens rotation matrix to zero out y
if y == 0
    c = 1;
    s = 0;
elseif abs(y) > abs(x)
    t = -x / y;
    s = 1 / sqrt(1 + t^2);
    c = s*t;
else
    t = -y / x;
    c = 1 / sqrt(1 + t^2);
    s = c*t;
end
end

该代码输入为一个$n\times n$的实对称矩阵$A$和一个收敛精度的阈值$tol$，输出为一个$n\times 1$的特征值向量和一个$n\times n$的特征向量矩阵。其中用到的函数 `hessenberg` 是将矩阵化为上 Hessenberg 形式的函数，`givens_rotation` 是计算 Givens 旋转矩阵的函数。在算法的实现中，采用了带 Wilkinson 位移的隐式对称 QR 算法，并在 QR 迭代中利用 Givens 旋转矩阵实现矩阵的相似变换。