Gauss消去法、列主元素消去法及LU分解法的MATLAB实现

去年的时候写的，整理一下发上来

Gauss消去法

代码：

function x=Gauss_elimination(A,b)
tic
[rows,~]=size(A);
aug_mat=[A,b];
for i=1:rows
    if aug_mat(i,i)~=0
        coefficient=aug_mat(:,i);
        aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
        coefficient=-coefficient;
        coefficient(i)=0;
        aug_mat=coefficient*aug_mat(i,:)+aug_mat;
    else
        aug_mat(:,rows+1:end)=nan;
        break
    end
end
x=aug_mat(:,rows+1:end);
toc
end

实例：
$$A=\left[\begin{array}{ccc}2& 2& 3\\ 4& 7& 7\\ -2& 4& 5\end{array}\right]b=\left[\begin{array}{c}3\\ 1\\ -7\end{array}\right]\text{\hspace{1em}(1)}$$
$$A=\left[\begin{array}{ccc}2& 2& 3\\ 4& 7& 7\\ -2& 4& 5\end{array}\right]b=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\text{\hspace{1em}(2)}$$
$$A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]b=\left[\begin{array}{c}1\\ 2\end{array}\right]\text{\hspace{1em}(3)}$$
$$A=\left[\begin{array}{cc}0.0001& 1\\ 1& 1\end{array}\right]b=\left[\begin{array}{c}1\\ 2\end{array}\right]\text{\hspace{1em}(3)}$$
实例代码：

A1=[2,2,3;4,7,7;-2,4,5];
b1=[3;1;-7];
x1=Gauss_elimination(A1,b1)

A2=[2,2,3;4,7,7;-2,4,5];
b2=eye(3);
x2=Gauss_elimination(A2,b2)

A3=[0,1;1,1];
b3=[1;2];
x3=Gauss_elimination(A3,b3)

A4=[0.0001,1;1,1];
b4=[1;2];
x4=Gauss_elimination(A4,b4)

实例运行结果：
历时 0.001997 秒。
$$x1=\left[\begin{array}{c}2\\ -2\\ 1\end{array}\right]$$
历时 0.002637 秒。
$$x2=\left[\begin{array}{ccc}0.1944& 0.0556& -0.1944\\ -0.9444& 0.4444& -0.0556\\ 0.8333& -0.3333& 0.1667\end{array}\right]$$
历时 0.000716 秒。
$$x3=\left[\begin{array}{c}NaN\\ NaN\end{array}\right]$$
历时 0.001793 秒。
$$x4=\left[\begin{array}{c}1.0001\\ 0.9999\end{array}\right]$$
我们发现gauss消去法并不能应对主对角线在过程中出现0的情况，如x3结果，因此我们在之后给出列主元素消去法及LU分解法代码。

列主元素消去法

代码：

function x=principal_element(A,b)
tic
[rows,~]=size(A);
aug_mat=[A,b];
for i=1:rows
    temp_list=aug_mat(:,i);
    temp_list(1:max(1,i-1))=0;
    [~,exchange_pos]=max(abs(temp_list));
    aug_mat([i,exchange_pos],:)=aug_mat([exchange_pos,i],:);
    coefficient=aug_mat(:,i);
    coefficient=-coefficient./coefficient(i);
    coefficient(i)=0;
    aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
divisor=aug_mat((1:rows)+(0:rows:(rows*(rows-1))))';
aug_mat=aug_mat./divisor;
x=aug_mat(:,rows+1:end);
toc
end

我们发现列主元消去法能够解决x3的情况：

A3=[0,1;1,1];
b3=[1;2];
x3=principal_element(A3,b3)

历时 0.020405 秒。
$$x3=\left[\begin{array}{c}1\\ 1\end{array}\right]$$

LU分解法

第一部分：矩阵LU分解

function [L,U]=LU(A)
[rows,~]=size(A);
temp_mat=A;
L=zeros(rows);
for i=1:rows
    coefficient=temp_mat(:,i);
    coefficient=coefficient./coefficient(i);
    coefficient(1:i)=0;
    L(:,i)=coefficient;
    temp_mat=-coefficient*temp_mat(i,:)+temp_mat;
end
U=temp_mat;
L(eye(rows)==1)=1;
end

第二部分：依据LU矩阵求解线性方程组

function x=LUsolve(L,U,b)
[rows,~]=size(L);
aug_mat=[L,b];
for i=1:rows
    aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
    coefficient=-aug_mat(:,i);
    coefficient(1:i)=0;
    aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
aug_mat=[U,aug_mat(:,rows+1:end)];
for i=rows:-1:1
    aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
    coefficient=-aug_mat(:,i);
    coefficient(i:end)=0;
    aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
x=aug_mat(:,rows+1:end);
end

使用方法：以问题(1)为例

A=[2,2,3;4,7,7;-2,4,5];
b=[3;1;-7];
[L,U]=LU(A)
x=LUsolve(L,U,b)

求解结果:
$$L=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -1& 2& 1\end{array}\right]U=\left[\begin{array}{ccc}2& 2& 3\\ 0& 3& 1\\ 0& 0& 6\end{array}\right]x=\left[\begin{array}{c}2\\ -2\\ 1\end{array}\right]$$