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

去年的时候写的，整理一下发上来Gauss消去法代码：function x=Gauss_elimination(A,b)tic[rows,~]=size(A);aug_mat=[A,b];for i=1:rowsif aug_mat(i,i)~=0coefficient=aug_mat(:,i);aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);coefficient=

slandarer

9851人浏览 · 2021-06-10 16:25:21

slandarer · 2021-06-10 16:25:21 发布

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

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=\begin{bmatrix} 2 & 2 & 3 \\ 4 & 7 & 7 \\ -2 & 4 & 5 \end{bmatrix} b=\begin{bmatrix} 3\\ 1 \\ -7 \end{bmatrix} \tag{1}$
$A=\begin{bmatrix} 2 & 2 & 3 \\ 4 & 7 & 7 \\ -2 & 4 & 5 \end{bmatrix} b=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \tag{2}$
$A=\begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} b=\begin{bmatrix} 1 \\ 2 \end{bmatrix} \tag{3}$
$A=\begin{bmatrix} 0.0001 & 1 \\ 1 & 1 \end{bmatrix} b=\begin{bmatrix} 1 \\ 2 \end{bmatrix} \tag{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=\begin{bmatrix} 2\\ -2\\ 1 \end{bmatrix}$
历时 0.002637 秒。
$x2=\begin{bmatrix} 0.1944 & 0.0556 & -0.1944\\ -0.9444 &0.4444 &-0.0556\\ 0.8333 &-0.3333&0.1667\\ \end{bmatrix}$
历时 0.000716 秒。
$x3=\begin{bmatrix} NaN\\ NaN \end{bmatrix}$
历时 0.001793 秒。
$x4=\begin{bmatrix} 1.0001\\ 0.9999 \end{bmatrix}$
我们发现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=\begin{bmatrix} 1\\ 1 \end{bmatrix}$

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=\begin{bmatrix} 1 &0 & 0\\ 2 &1 &0\\ -1 &2 & 1\\ \end{bmatrix} U=\begin{bmatrix} 2 & 2 & 3\\ 0 & 3 & 1\\ 0 & 0 & 6\\ \end{bmatrix} x=\begin{bmatrix} 2\\ -2\\ 1 \end{bmatrix}$