% 幂法计算矩阵最大特征值及对应特征向量（脚本版）
% 功能：无需函数调用，直接输入矩阵A和初始向量u，输出迭代过程、收敛后的最大特征值及特征向量
% 运行方式：直接运行此m脚本，无需命令行操作，结果自动输出
clear; clc;  % 清空工作区变量与命令行窗口，避免历史数据干扰


%% 1. 输入待计算的矩阵A与初始特征向量u
% 示例1：3阶实矩阵（可替换为任意方阵，确保存在唯一最大特征值）
A = [4 1 1; 
     1 3 2; 
     1 2 3];  
% 示例2：对称矩阵（幂法对对称矩阵收敛更稳定，可替换上方A）
% A = [17 3 4; 3 1 12; 4 12 8];

% 初始特征向量u（需为非零向量，通常设为全1列向量，维度与A匹配）
u = [1; 1; 1];  % 注意：需为列向量（满足矩阵乘法A*u的维度规则）

fprintf('待计算最大特征值的矩阵A：\n');
disp(A);
fprintf('\n初始特征向量u：\n');
disp(u);
[m,n] = size(A);  % 获取矩阵A的维度（需满足m=n，即方阵）


%% 2. 幂法核心计算逻辑（整合原Matrix_pow函数核心代码）
i = 1;                % 迭代次数计数器（初始为1）
dellab = 1;           % 相邻两次特征值的差值（初始设为1，满足循环条件）
miu(1) = 1;           % 特征值迭代序列（初始值设为1）
u_normal = [];        % 存储迭代过程：[迭代次数, 特征向量元素]

% 迭代循环：当特征值差值小于1e-6时收敛
while abs(dellab) > 1e-6
    % 记录当前迭代次数与特征向量（第1列是迭代次数，后续列是特征向量元素）
    u_normal = [u_normal; i-1, u'];  % i-1为迭代次数（从0开始计数）
    % 步骤1：计算v = A*u（特征向量迭代）
    v = A * u;
    % 步骤2：找v中绝对值最大的元素，确定当前特征值miu(i)
    [vmax, j] = max(abs(v));        % vmax是v的最大绝对值，j是对应索引
    miu(i) = vmax * sign(v(j));     % 特征值（保留符号，确保正负正确）
    % 步骤3：归一化特征向量u（避免数值溢出）
    u = v / miu(i);
    % 步骤4：计算相邻特征值的差值，判断是否收敛
    if i >= 2  % 从第2次迭代开始计算差值
        dellab = miu(i) - miu(i-1);
    end
    i = i + 1;  % 迭代次数加1
end

% 最终结果赋值（与原函数输出对应）
lambta1 = miu;                  % 所有迭代步骤的特征值序列
final_lambta = lambta1(end);    % 收敛后的最大特征值
final_u = u;                    % 收敛后的对应特征向量


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

fprintf('\n2. 特征值迭代序列（最后一个值为收敛后的最大特征值）：\n');
disp(lambta1);

fprintf('\n3. 收敛后的最大特征值 lambta1：%.8f\n', final_lambta);
fprintf('\n4. 对应最大特征值的特征向量 u（归一化后）：\n');
disp(final_u);

fprintf('\n5. 迭代收敛条件：相邻特征值差值 |dellab| ≤ 1e-6\n');
fprintf('   最终迭代次数：%d 次\n', length(lambta1)-1);  % 迭代次数=特征值序列长度-1


%% 4. 结果验证（验证 A*u ≈ lambta1*u，确认特征值/向量正确性）
fprintf('\n==================== 结果验证 ====================\n');
left = A * final_u;              % 计算 A*u
right = final_lambta * 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('说明：误差为浮点计算误差，数值极小时表明特征值/向量正确。\n');