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');