matlab 主因素高斯消除 全组元消去法的实现

在确定精度下（即只保留一定小数计算的情况下），在解多元方程组时，可能会遇到大数吃小数导致误差放大的情况。因此，为了避免出现大数吃小数的情况，在求解过程中需要将主元（大数）移至前端进行处理。

可以参考《数值计算方法》这本书中的其中例题，加深理解

 

%高斯消去法+主因素
%author LijiaYi（foddcus） FaFu  2022.3
%email：foddcus@163.com
%最后修改日期2022-5-1


%function Toutput=MAeliminationM(input);%模块化功能

%input输入矩阵，包括系数和参数
%oringinData保留的初始矩阵
%output 输出答案，对应X的各个值
%Toutput：修正顺序后的答案，对应输入数据的X顺序
%rememberdet：行列变化中介值
%remember：行列变化记录矩阵
%xnum,ynum 输入矩阵的尺寸
%Lneed 消元中，主式需要乘的变量
%keepnum 保留小数位数
%contrastFrom（各个多项式的误差）
%simdata 模拟序列 用来修正输出序列
%deltanum：矩阵中变化的最大值

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

input=[1,3,4,3;
    4,1,3,1;
    4,4,6,11];%输入矩阵
oringinData=input;



error1='the condition are not meet for using this algorithm'

%format short%设置精度 short为精确到小数点后3位 数值型 roundn（A，3）保留3位小数

keepnum=-10;%设置保留n位小数,记得加负号表示小数位
[ynum,xnum]=size(input);%获取输入矩阵大小

%创建模拟序列,用来修正顺寻
for i=1:ynum
    simData(i,1)=i
end

if ynum==xnum-1; %满足高斯消去法的要求
remember=[];
    for i=1:ynum-1
        %         for k=i:ynum %这里用朴素的方法找到剩余矩阵的最大值,当然也可以使用max函数求，不过这就相当于用工具箱了
        %             for u=i:ynum
        deltanum=input(i,i)%%先做主因素提取
        %                 [bigNum,L]=max(input);
        %                 [bigNum,Lbing]=max(bigNum);
        %                 realR=L(Lbing);
        %                 realL=Lbing;
        %             end
        %         end
        
        for k=i:ynum %这里用朴素的方法找到剩余矩阵的最大值,当然也可以使用max函数求，不过这就相当于用工具箱了
            for u=i:ynum
                if deltanum<input(k,u)
                    deltanum=input(k,u);
                    remember(i,1:3)=[k,u,i] ;%算出最大位置所在
                end
            end
        end
        [CEK1,~]=size(remember);%如果没有找到更大的值，使用原位置
        if CEK1<i
            remember(i,1:3)=[i,i,i];
        end
        remember(i,1:3)=[i,i,i];
        %这里开始做行列式变化
        k=remember(i,1);
        u=remember(i,2);
        m=remember(i,3)
        rememberdet=input(k,:);%行变化
        input(k,:)=input(i,:);
        input(i,:)=rememberdet;
        rememberdet=input(:,u);%列变化
        input(:,u)=input(:,i);
        input(:,i)=rememberdet;
        rememberdet=oringinData(:,u);%原始数据也来列变化，方便验算
        oringinData(:,u)=oringinData(:,i);
        oringinData(:,i)=rememberdet;
        rememberdet=simData(u,1);%模拟数据也来列变化，方便返还结果
        simData(u,1)=simData(i,1);
        simData(i,1)=rememberdet;
        
        
        for j=i+1:ynum
            if input(i,i)==0%检查是否符合要求
                disp(error1);
                break
            end
            
            Lneed(j,i)=input(j,i)/input(i,i);
            Lneed(j,i)=roundn(Lneed(j,i),keepnum);%四舍五入，下同
        end%做完系数的工作
        
        
        for j=i+1:ynum
            input(j,:)=input(j,:)-input(i,:).*Lneed(j,i);
            
        end%消元
    end
    
    
    %回代
    for i=1:ynum%这里计算出结果
        output(ynum+1-i,1)=roundn(input(ynum+1-i,xnum),keepnum)/roundn(input(ynum+1-i,xnum-i),keepnum);%output为输出结果
        output(ynum+1-i,1)=roundn(output(ynum+1-i,1),keepnum);
        for j=i+1:ynum
            %回代结果到多项式中，将其他值消除
            input(ynum+1-j,xnum)=roundn(input(ynum+1-j,xnum),keepnum)-roundn(input(ynum+1-j,xnum-i),keepnum)*output(ynum+1-i,1);
            input(ynum+1-j,xnum)=roundn(input(ynum+1-j,xnum),keepnum);
            input(ynum+1-j,xnum-i)=0;
        end
        
    end
    
    %验算
    contrast=0;
    for i=1:ynum
        for j=1:ynum
            contrast=contrast+oringinData(i,j)*output(j,1);
        end
        contrastForm(i,1)=roundn(contrast,keepnum);
        contrast=oringinData(i,end)-contrast;
        contrastForm(i,2)=roundn(contrast,keepnum);
        contrast=0;
    end
    
    for i=1:ynum%修正序列
        Toutput(simData(i,1),1)=output(i,1);
    end
    
else
    disp(error1);
end