clear; clc; 

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

% 原点平移参数p（作用：加速收敛，通常根据A的对角线元素估计，如取对角线均值/某一元素）
% 示例：A对角线元素为5、6、7，取p=6（接近特征值分布中心，可减少迭代次数）
p = 6;  

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

% 输出输入信息，便于对比
fprintf('待计算最大特征值的矩阵A：\n');
disp(A);
fprintf('\n原点平移参数 p = %.2f\n', p);
fprintf('\n初始特征向量u：\n');
disp(u);
[m,n] = size(A);  % 获取矩阵A维度（需满足m=n，即方阵）


%% 2. 原点平移幂法核心计算逻辑（整合原Matrix_parallel_trans_pow函数核心代码）
I = eye(n);               % 构造n阶单位矩阵
B = A - p * I;            % 原点平移后的矩阵：B = A - p*I

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

% 迭代循环：当特征值差值小于1e-6时收敛
while abs(dellab) > 1e-6
    % 记录当前迭代次数与特征向量（第1列=迭代次数，后续列=特征向量元素，迭代次数从0开始）
    u_normal = [u_normal; i-1, u'];  % u'将列向量转为行向量，便于拼接
    % 步骤1：特征向量迭代计算：v = B*u
    v = B * u;
    % 步骤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 >= 2
        dellab = miu(i) - miu(i-1);
    end
    i = i + 1;  % 迭代次数加1
end

% 最终结果计算（与原函数输出对应）
lambta1 = miu + p;                % A的特征值序列（B的特征值+平移参数p）
final_lambta_B = miu(end);        % B收敛后的最大特征值
final_lambta_A = lambta1(end);    % A收敛后的最大特征值（最终结果）
final_u = u;                      % 收敛后的对应特征向量（A和B的特征向量相同）


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

% 2. 特征值迭代序列对比（B的特征值 vs A的特征值）
fprintf('\n2. 特征值迭代序列对比：\n');
fprintf('   B = A - p*I 的特征值序列：\n');
disp(miu);
fprintf('   A 的特征值序列（B的特征值 + p）：\n');
disp(lambta1);

% 3. 最终核心结果
fprintf('\n3. 收敛后的核心结果：\n');
fprintf('   - B = A - p*I 的最大特征值：%.8f\n', final_lambta_B);
fprintf('   - A 的最大特征值 lambta1：%.8f\n', final_lambta_A);
fprintf('   - 对应最大特征值的特征向量 u（归一化后）：\n');
disp(final_u);

% 4. 迭代收敛信息
fprintf('\n4. 迭代收敛信息：\n');
fprintf('   - 收敛条件：相邻特征值差值 |dellab| ≤ 1e-6\n');
fprintf('   - 总迭代次数：%d 次\n', length(miu)-1);  % 迭代次数=特征值序列长度-1
fprintf('   - 原点平移作用：通过调整p使B的特征值分布更集中，加速收敛（对比普通幂法迭代次数更少）\n');


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

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