matlab 使用基本函数还原反幂法迭代 求相近的特征根

%&subnorm
%反幂法 
%author（foddcus） FaFu  2022.3 
%email：foddcus@163.com
%%
%输入：目标矩阵（input）；目标接近的特征值（charnum）、迭代上限（aimnum）、目标误差（errorS）
%输出：目标特征值（output）；xform、zfrom，包含x、z特征向量的值。最接近的特征值

%%
%function Tresult=eliminationM(input);%模块化功能

%input输入矩阵，包括系数和参数
%inputC 用于分解LU矩阵的合成矩阵
%oringinData保留的初始矩阵
%output 输出答案，对应X的各个值
%Toutput：修正顺序后的答案，对应输入数据的X顺序
%rememberdet：行列变化中介值
%remember：行列变化记录矩阵
%xnum,ynum 输入矩阵的尺寸
%Lneed 消元中，主式需要乘的变量
%keepnum 保留小数位数
%contrastFrom（各个多项式的误差）
%deltanum：矩阵中变化的最大值
%Ldet：分解出来的下矩阵
%Udet:分解出来的上矩阵
%Zform：z特征向量值表格
%xform：包含x特征向量的表格

clear all%启用模块化后需要注释该句

input=[1,2,3;
       1,0,2
       0,1,-1];%输入矩阵
charNum=-3; %charnum范围一般比最大特征根要大，比最小特征根要小，但差距不能过大 ，否则需要提高迭代上限来弥补精度
aimnum=400;
errorS=10^(-5);
%%
oringinData=input;
error1='the condition are not meet for using this algorithm' ;
%format short%设置精度 short为精确到小数点后3位 数值型 roundn（A，3）保留3位小数
keepnum=-20;%设置保留n位小数,记得加负号表示小数位
[ynum,xnum]=size(input);%获取输入矩阵大小
if ynum~=xnum
    disp(error1)
end

Idet=charNum*eye(ynum);
input=input-Idet;

for i=1:ynum
    for j=1:ynum
        Udet(i,j)=0;
        Ldet(i,j)=0;
    end
end
inputC=input


for i=1:ynum%构造合成矩阵（行j，列I）
    for j=i:ynum  %行逻辑为主，列逻辑为辅
        for I=i:xnum%I这里做判断的算子
%             if I>j%这里需要提前做一个判断，用于加减的转变
%                 overnum=j;
%             else
%                 overnum=I;
%             end
     
            if j==i%跳过该行
            else
                if I==i
                    inputC(j,I)=roundn(inputC(j,I)/inputC(i,i),keepnum);
                else
                    inputC(j,I)=roundn(inputC(j,I)-inputC(i,I)*inputC(j,i),keepnum);
                end   
            end         
        end
    end    
end

for i=1:xnum%将合成矩阵拆分为L、U矩阵
    for j=1:ynum
        if j==i
            Ldet(j,i)=1;
            Udet(j,i)=inputC(j,i);
        end
        if j>i
            Ldet(j,i)=inputC(j,i);
        end
        if j<i&&i<4
            Udet(j,i)=inputC(j,i);
        end
    end
end

%%开始计算

X=[0;0;1];
Y=[0;0;1];       
alpa(1,1)=1;
beta(1,1)=1
for k=1:aimnum
    Y=X/alpa(k,:);
    Z=(Ldet^(-1))*Y;
    X=(Udet^(-1))*Z;
    zform(k,:)=Z;
    xform(k,:)=X;
    alpa(k+1,1)=max(max(X));
    fiv(k,1)=charNum+1/alpa(k+1,1);
   if errorS>abs(1/alpa(k+1,1)-1/alpa(k,1))
         break
   end
end
disp(["临近特征值",num2str(fiv(k,1))]);