数值分析·学习 | 解线性方程组的直接方法（高斯消去法以及LU求解）matlab实现

一、前言：

个人学习内容分享

二、算法描述：

1、高斯消去法：

设有线性方程组

$\left\{\begin{matrix} a_{11}x_{1}+a_{12}x_{2}+\cdots +a_{1n}x_{n}=b_{1},\\ a_{21}x_{1}+a_{22}x_{2}+\cdots +a_{2n}x_{n}=b_{2},\\ \vdots \\ a_{n1}x_{1}+a_{n2}x_{2}+\cdots +a_{nn}x_{n}=b_{n}, \end{matrix}\right.$ 编辑

或写为矩阵形式

$\left(\begin{array}{cccc} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1 n}^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2 n}^{(2)} \\ & & \ddots & \vdots \\ & & & a_{n n}^{(n)} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right)=\left(\begin{array}{c} b_{1}^{(1)} \\ b_{2}^{(2)} \\ \vdots \\ b_{n}^{(n)} \end{array}\right)$ 编辑

如果 $a_{kk}^{(k)}\neq 0(k=1,2,\cdots ,n)$ 编辑，则可通过高斯消去法将Ax=b约化为等价的三角形线性方程组，形式如下：

计算公式如下：

消元计算(k=1,2,···,n-1)

$\left\{\begin{matrix} m_{ik}=a_{ik}^{(k)}/a_{kk}^{(k)}, & i=k+1,\cdots,n,\\ a_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ik}a_{kj}^{(k)},&i,j=k+1,\cdots,n,\\ b_{i}^{(k+1)}=b_{i}^{(k)}-m_{ik}b_{k}^{(k)},&i=k+1,\cdots,n, \end{matrix}\right.$ 编辑

回代计算

$\left\{\begin{matrix} x_{n}=b_{n}^{(n)}/a_{nn}^{(n)},\\ x_{i}=(b_{i}^{(i)}-\sum_{j=i+1}^{n}a_{ij}^{(i)}x_{j})/a_{ii}^{(i)},i=n-1,\cdots ,2,1 \end{matrix}\right.$ 编辑

2、高斯消去法-列主元消去法：

1.det←1

2.对于k=1,2,···,n-1

（1）按列选主元

$\left | a_{i_{k},k} \right |=\underset{k\leqslant i \leqslant n}{max}\left | a_{ik} \right |$ 编辑

（2）如果 $a_{i_{k},k}=0$ 编辑,则停止计算(det(A)=0)

（3）如果 $i_{k}=k$ 编辑则转（4）

换行： $a_{kj}\leftrightarrow a_{i_{k},j}(j=k,k+1,\cdots ,n)$ 编辑

$b_{k}\leftrightarrow b_{i_{k}}$ 编辑

$det\leftarrow -det$ 编辑

（4）消元计算

对于i=k+1,···,n

$a_{ik}\leftarrow m_{ik}=a_{ik}/a_{kk}$ 编辑

对于j=k+1,···,n

$a_{ij}\leftarrow a_{ij}-m_{ik}*a_{kj}$ 编辑

$b_{i}\leftarrow b_{i}-m_{ik}/b_{k}$ 编辑

（5） $det\leftarrow a_{kk}*det$ 编辑

3.如果 $a_{nn}=0$ 编辑，则计算停止（det(A)=0）

4.回代求解

（1） $b_{n}\leftarrow b_{n}/a_{nn}$ 编辑

（2）对于i=n-1,···,2,1

$b_{i}\leftarrow (b_{i}-\sum_{j=i+1}^{n}a_{ij}*b_{j})/a_{ii}$ 编辑

5. $det\leftarrow a_{nn}*det$ 编辑

3、LU分解：

$u_{1i}=a_{1i},(i=1,2,\cdots ,n),l_{i1}=\frac{a_{i1}}{u_{11}},i=2,3,\cdots ,n.$ 编辑

计算U的第r行，L的第r列元素（r=2,3,...,n）:

$u_{ri}=a_{ri}-\sum_{k=1}^{r-1}l_{rk}y_{k},i=r,r+1,\cdots ,n;$ 编辑

$L_{ir}=\frac{(a_{ir}-\sum_{k=1}^{r-1}l_{ik}u_{kr})}{u_{rr}},i=r+1,\cdots ,n,and\: r\neq n$ 编辑

求解Ly=b,Ux=y的计算公式：

$\begin{cases} & \text{ }y_{1}=b_{1} \\ & \text{ }y_{i}=b_{i}-\sum_{k=1}^{i-1}l_{ik}y_{k} ,i=2,3,\cdots ,n; \end{cases}$ 编辑

$\begin{cases} & \text{ }x_{n}=\frac{y_{n}}{u_{nn} } \\ & \text{ }x_{i}=\frac{y_{i}-\sum_{k=i+1}^{n}u_{ik}x_{k}}{u_{ii}} ,i=n-1,n-2,\cdots ,1; \end{cases}$ 编辑

三、实现代码：

1、高斯消去法：

function x=Gauss_Elimination_Method(A)
%功能：一般的高斯消元法
%依次消元
for k=1:length(A)-2
    for m=1:length(A)-1-k
        A(k+m,:)=A(k+m,:)-A(k+m,k)/A(k,k).*A(k,:);
    end
end
%回代求解X
x=zeros(length(A)-1,1);
for k=1:length(A)-1
    temp=length(A)-k;
   x(temp,1)=(A(temp,length(A))-sum(A(temp,1:length(A)-1)*x(:,1)))/A(temp,temp);
end
end

2、高斯消去法-列主元消去法：

function x=Gauss_Elimination_Method_GAI(A)
%功能：列主元的高斯消元法
%依次消元
for k=1:length(A)-2
       temp1=max(A(k:length(A)-1,k));
     for m=k:length(A)-1
         if A(m,k)==temp1
             temp1=A(m,:);
             A(m,:)=A(k,:);
             A(k,:)=temp1;
             continue;
         end
     end
    for m=1:length(A)-1-k
        A(k+m,:)=A(k+m,:)-A(k+m,k)/A(k,k).*A(k,:);
    end
end
%回代求解X
x=zeros(length(A)-1,1);
for k=1:length(A)-1
    temp3=length(A)-k;
   x(temp3,1)=(A(temp3,length(A))-sum(A(temp3,1:length(A)-1)*x(:,1)))/A(temp3,temp3);
end
end

3、LU分解：

function x=LU(A,b)
%功能：LU算法
L=eye(length(A));
for k=1:length(A)-1
    for m=1:length(A)-k
        L(k+m,k)=A(k+m,k)/A(k,k);
        A(k+m,:)=A(k+m,:)-A(k+m,k)/A(k,k).*A(k,:);
    end
end
y=L\b;
x=A\y;
end

4、求逆矩阵：

function X=Inverse_matrix(A)
%功能：求逆矩阵
L=length(A);
I=eye(L);
A=[A,I];
for k=1:L-1
    for m=1:L-k
        A(k+m,:)=A(k+m,:)-A(k+m,k)/A(k,k).*A(k,:);
    end
end
for k=1:L-1
    for m=L-k:-1:1
        if A(m+1,:)==0
            continue;
        else
            A(m,:)=A(m,:)-A(m,L+1-k)/A(L+1-k,L+1-k).*A(L+1-k,:);
        end
    end
end
for k=1:L
    A(k,:)=A(k,:)/A(k,k);
end
X=A(:,L+1:2*L);
end