Jacobi.m

% 优化版Jacobi方法（解决卡住问题）
% 新增：最大迭代次数限制、非对称矩阵强制退出、索引有效性检查
clear; clc;

%% 1. 输入参数（可修改）
A = [4 1 1; 
     1 3 2; 
     1 2 3];  % 对称矩阵（示例）
% A = [1 2 3; 4 5 6; 7 8 9];  % 非对称矩阵（测试用，会快速退出）
X = 1;  % 1=上三角模式，2=下三角模式
max_iter = 10000;  % 最大迭代次数（防止无限循环）
tol = 1e-8;        % 收敛阈值


%% 2. 初始化与参数检查
[m,n] = size(A);
if m ~= n
    error('矩阵必须为方阵！');
end
% 检查对称性，非对称矩阵提示并限制迭代
is_symmetric = isequal(A, A');
if ~is_symmetric
    warning('输入矩阵非对称，Jacobi方法可能不收敛，将限制迭代次数！');
end

C2 = 1;          % 最大非对角元素绝对值
Rk = eye(n);     % 旋转矩阵累积器
k = 0;           % 迭代计数器

% 选择旋转模式（上三角/下三角）
switch X
    case 1
        fun = @(M) triu(M,1);  % 提取上三角非对角元素（不含对角线）
    case 2
        fun = @(M) tril(M,-1); % 提取下三角非对角元素（不含对角线）
    otherwise
        error('X必须为1或2！');
end


%% 3. 核心迭代（添加退出保护）
while abs(C2) > tol && k < max_iter
    % 提取非对角元素并找最大值位置
    A1 = fun(A);
    [C1, I1] = max(abs(A1));  % 每列最大绝对值及行索引
    [C2, j] = max(C1);        % 全局最大绝对值及列索引
    i = I1(j);                % 行索引
    
    % 保护：避免i=j（对角线元素）或索引越界
    if i == j || i < 1 || i > n || j < 1 || j > n
        warning('无效索引(i,j)，提前退出迭代！');
        break;
    end
    
    % 计算Givens旋转参数（避免除以零）
    a_ij = A(i,j);
    if abs(a_ij) < 1e-15  % 防止分母为零
        C2 = 0;  % 视为已收敛
        continue;
    end
    dij = (A(i,i) - A(j,j)) / (2 * a_ij);
    t = sign(dij) / (abs(dij) + sqrt(1 + dij^2));
    c = 1 / sqrt(1 + t^2);
    s = t * c;
    
    % 构造旋转矩阵并更新A和Rk
    R = eye(n);
    R(i,i) = c;  R(i,j) = s;
    R(j,i) = -s; R(j,j) = c;
    A = R * A * R';  % 相似变换
    Rk = R * Rk;     % 累积旋转
    k = k + 1;       % 计数迭代次数
end


%% 4. 输出结果与状态判断
fprintf('迭代结束：\n');
if k >= max_iter && abs(C2) > tol
    fprintf('警告：达到最大迭代次数（%d次），未完全收敛！\n', max_iter);
else
    fprintf('成功收敛：迭代次数=%d次，最终非对角最大元素=%.2e\n', k, C2);
end

fprintf('\n对角化矩阵Ak：\n');
disp(A);
fprintf('\n正交矩阵Qk：\n');
disp(Rk');


%% 5. 结果验证
fprintf('\n验证：\n');
fprintf('Qk的正交性（Qk''*Qk）：\n');
disp(Rk' * Rk);
fprintf('原始矩阵恢复（Qk*Ak*Qk''）：\n');
disp(Rk * A * Rk');