clear; clc; 


%% 1. 输入计算所需参数（可根据需求替换）
% 待求特征值的方阵A（示例为3阶对称矩阵，确保B=A-pI可逆，可替换为任意方阵）
A = [6 2 1; 
     2 5 3; 
     1 3 7];  

% 原点平移参数p（目标：求解A在p附近的特征值，若求最小特征值，可设p接近0或根据特征值估计值设置）
% 示例：假设A的特征值在1~10之间，取p=2（靠近较小特征值，加速收敛到最小特征值）
p = 2;  

% 初始特征向量u（需为非零列向量，维度与A匹配，通常设为全1列向量）
u = [1; 1; 1];  

% 输出输入信息，便于后续对比
fprintf('待计算特征值的矩阵A：\n');
disp(A);
fprintf('\n原点平移参数 p = %.2f（求解A在p附近的特征值）\n', p);
fprintf('\n初始特征向量u：\n');
disp(u);
[m,n] = size(A);  % 获取矩阵A维度（需满足m=n，即方阵）
if m ~= n
    error('输入矩阵A必须为方阵！');  % 反幂法仅适用于方阵
end


%% 2. 反幂法（带原点平移）核心计算逻辑（整合原Matrix_invpow2函数核心代码）
I = eye(n);               % 构造n阶单位矩阵
B = A - p * I;            % 原点平移后的矩阵：B = A - pI（核心：通过B⁻¹的特征值反推A的特征值）

i = 1;                    % 迭代次数计数器（初始为1）
dellab = 1;               % 相邻特征值差值（初始设为1，满足循环条件）
miu(1) = max(u);          % B⁻¹的特征值迭代序列（初始值设为u的最大元素）
u_normal = [];            % 存储迭代过程：[迭代次数, 特征向量元素]

% 迭代循环：当特征值差值小于1e-6时收敛
while abs(dellab) > 1e-6
    % 记录当前迭代次数与特征向量（第1列=迭代次数，后续列=特征向量元素，迭代次数从0开始）
    u_normal = [u_normal; i-1, u'];  % u'将列向量转为行向量，便于拼接
    % 步骤1：反幂法核心（解线性方程组Bv=u，等价于v = B⁻¹u，避免直接求逆矩阵）
    v = B \ u;  % 利用MATLAB的线性方程组求解器，比直接求逆更稳定
    i = i + 1;  % 迭代次数加1（因v是下一步结果，此处先更新i）
    % 步骤2：确定B⁻¹当前迭代的特征值miu(i)（取v中绝对值最大元素，保留符号）
    [vmax, j] = max(abs(v));        % vmax=v的最大绝对值，j=对应索引
    miu(i) = vmax * sign(v(j));     % B⁻¹的特征值（符号由v(j)决定）
    % 步骤3：归一化特征向量u（避免数值溢出，确保迭代稳定）
    u = v / miu(i);
    % 步骤4：计算相邻特征值差值（从第2次迭代开始判断收敛）
    if i >= 3  % 因i已在v计算后加1，此处i=3对应第2次迭代的差值
        dellab = miu(i) - miu(i-1);
    end
end

% 最终结果计算（与原函数输出对应）
lambta = 1 ./ miu + p;             % A的特征值序列（λ = 1/μ + p，μ为B⁻¹的特征值）
final_miu = miu(end);              % B⁻¹收敛后的特征值
final_lambta = lambta(end);        % A收敛后的目标特征值（最终结果）
final_u = u;                       % 收敛后的对应特征向量（A和B⁻¹的特征向量相同）


%% 3. 输出计算结果
fprintf('\n==================== 反幂法（带原点平移）计算结果 ====================\n');
% 1. 迭代过程记录
fprintf('1. 迭代过程（第1列：迭代次数，后续列：特征向量）：\n');
disp(u_normal);

% 2. 特征值迭代序列关系（B⁻¹的特征值μ vs A的特征值λ=1/μ+p）
fprintf('\n2. 特征值迭代序列关系：\n');
fprintf('   B⁻¹ = (A - pI)⁻¹ 的特征值序列μ：\n');
disp(miu);
fprintf('   A 的特征值序列λ = 1/μ + p：\n');
disp(lambta);

% 3. 最终核心结果
fprintf('\n3. 收敛后的核心结果：\n');
fprintf('   - B⁻¹ 的特征值μ：%.8f\n', final_miu);
fprintf('   - A 在p附近的特征值λ：%.8f\n', final_lambta);
fprintf('   - 对应特征值的特征向量u（归一化后）：\n');
disp(final_u);

% 4. 迭代收敛信息
fprintf('\n4. 迭代收敛信息：\n');
fprintf('   - 收敛条件：相邻特征值差值 |dellab| ≤ 1e-6\n');
fprintf('   - 总迭代次数：%d 次\n', length(miu)-1);  % 迭代次数=特征值序列长度-1
fprintf('   - 反幂法优势：比普通幂法更快收敛到指定p附近的特征值（尤其适合求最小特征值）\n');


%% 5. 结果验证（验证 A*u ≈ λu，确认特征值/向量正确性）
fprintf('\n==================== 结果验证 ====================\n');
left = A * final_u;              % 计算 A*u（特征方程左侧）
right = final_lambta * final_u;  % 计算 λ*u（特征方程右侧）
diff = left - right;             % 两者差值（理论应为0向量，误差源于浮点计算）
error = norm(diff, 'fro');       % 用Frobenius范数量化误差（越小越准确）

fprintf('A*u 的计算结果：\n');
disp(left);
fprintf('λ*u 的计算结果：\n');
disp(right);
fprintf('\n验证误差（||A*u - λ*u||_F）：%.6e\n', error);
fprintf('说明：误差为浮点计算误差，数值≤1e-6时表明特征值/向量正确。\n');