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