高斯消元法和列主元素法

高斯消元法

设线性方程组$Ax=b$，其中
$$(A^{(1)},b^{(1)})=\left(\begin{array}{cccccc}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_1^{(1)}\\ a_{21}^{(1)}& a_{22}^{(1)}& a_{23}^{(1)}& \dots & a_{2n}^{(1)}& b_2^{(1)}\\ a_{31}^{(1)}& a_{32}^{(1)}& a_{33}^{(1)}& \dots & a_{3n}^{(1)}& b_3^{(1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}^{(1)}& a_{n2}^{(1)}& a_{n3}^{(1)}& \dots & a_{nn}^{(1)}& b_n^{(1)}\end{array}\right)$$
用高斯消元法时，步骤如下：

1. 计算$-l_{i1}=-\frac{a_{i1}^{(1)}}{a_{11}^{(1)}},(i=2,3,...,n)$，然后将矩阵的第一行乘这个系数分别加到每一行，这样，第一列的对角元一下就全是0。
   $$\left(A^{(2)},b^{(2)}\right)=\left(\begin{array}{cccccc}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_1^{(1)}\\ 0& a_{22}^{(2)}& a_{23}^{(2)}& \dots & a_{2n}^{(2)}& b_2^{(2)}\\ 0& a_{32}^{(2)}& a_{33}^{(2)}& \dots & a_{3n}^{(2)}& b_3^{(2)}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& a_{n2}^{(2)}& a_{n3}^{(2)}& \dots & a_{nn}^{(2)}& b_n^{(2)}\end{array}\right]$$
   然后用除去第一行和第一列的矩阵重复上述步骤，n-1步之后，将会得到
   $$\left(A^{(n)},b^{(n)}\right)=\left(\begin{array}{llllll}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_1^{(1)}\\ & a_{22}^{(2)}& a_{23}^{(2)}& \dots & a_{2n}^{(2)}& b_2^{(2)}\\ & & a_{33}^{(3)}& \dots & a_{3n}^{(3)}& b_3^{(3)}\\ & & & \ddots & ⋮& ⋮\\ & & & & a_{nn}^{(n)}& b_n^{(n)}\end{array}\right)$$

这个时候消元完成，然后进行求解，如下：
$$\{\begin{array}{l}{x}_n=b_n^{(n)}/a_{nn}^{(n)}\\ x_i=\left(b_i^{(i)}-{\sum }_{j=i+1}^n{a}_{ij}^{(i)}{x}_j\right)/a_{ii}^{(i)},i=n-1,n-2,\cdots ,1\end{array}$$

matlab代码如下：

function x=gaosi(A,b)
[rows,cols]=size(A);
if rows~=cols
    error("rows != cols")
end
n=rows;
Ab=[A,b];
%消元
for i=1:n-1
   l= Ab((i+1):end,i)/Ab(i,i);
   rep_Ab_row = repmat(Ab(i,:),n-i,1);
   Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;
end
A=Ab(:,1:n);
b=Ab(:,n+1);
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
   x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);
end

关于代码中的l.*rep_Ab_row，l是一个m行1列的系数矩阵，rep_Ab_row是一个m行n+1列的矩阵，这里有点像numpy里的广播机制，前者每一行乘到后者每一行的每个元素上。
比如

>> a

a =

     1
     2
     3

>> b

b =

     1     2
     3     4
     5     6

>> a.*b

ans =

     1     2
     6     8
    15    18

列主元素法

高斯消元法虽然简单，但是要确保消元过程中$a_{i,i}^{(i)}\ne 0$，因为要作为除数，而且当它比较小的时候，作为除数会导致结果特别大，使得计算过程中的误差变大。
列主元素法就是在每一次消元之前，先对A进行初等行变换，即在这一列中，选取绝对值最大元素所在的行，然后和当前的第一行交换，可以使得误差较小。

%列主元素法
function x=lzys(A,b)
[rows,cols]=size(A);
if rows~=cols
    error("rows != cols")
end
n=rows;
Ab=[A,b];
%消元
for i=1:n-1
   %先找列主元素
   [M,idx] = max(abs(Ab(i:end,i)));
   %交换
   idx = idx + i-1;%max得到的idx=1时，指的是第i行
   Ab([i,idx],:)=Ab([idx,i],:);%去掉这个分号，可以选择列主元素的结果
   l= Ab((i+1):end,i)/Ab(i,i);
   rep_Ab_row = repmat(Ab(i,:),n-i,1);
   Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;
end
A=Ab(:,1:n);
b=Ab(:,n+1);
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
   x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);
end

改动

在求不同行的系数乘第一行时，之前采用的是一个笨方法，即先把当前第一行重复复制为(n-i)行，然后对每行乘那个对应的系数。

rep_Ab_row = repmat(Ab(i,:),n-i,1);
Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;

之后发现这个过程可以直接由矩阵相乘得到，计算的系数矩阵是一个m行1列的矩阵，可以直接和1行n+1列的矩阵相乘，效果和前边先重复，再单独乘一样。

Ab((i+1):end,:) = Ab((i+1):end,:)-l*Ab(i,:)

同时，循环中，只需要对右下角的矩阵进行处理，我之前的代码中对所有列都进行处理了，虽然当前矩阵左边的列都是0。这个也可以改动
改动后的代码

%高斯消元法
function x=gaosi(A,b)
[rows,cols]=size(A);
if rows~=cols
    error("rows != cols")
end
n=rows;
Ab=[A,b];
%消元
for i=1:n-1
   l= Ab((i+1):end,i)/Ab(i,i);
   %rep_Ab_row = repmat(Ab(i,:),n-i,1);
   %上边两步可以用矩阵乘法替代
   Ab((i+1):end,i:end) = Ab((i+1):end,i:end)-l*Ab(i,i:end);
   %Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;
end
A=Ab(:,1:n);
b=Ab(:,n+1);
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
   x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);
end

%列主元素法
function x=lzys(A,b)
[rows,cols]=size(A);
if rows~=cols
    error("rows != cols")
end
n=rows;
Ab=[A,b];
%消元
for i=1:n-1
   %先找列主元素
   [M,idx] = max(abs(Ab(i:end,i)))
   %交换
   idx = idx + i-1;%max得到的idx=1时，指的是第i行
   Ab([i,idx],i:end)=Ab([idx,i],i:end);
   l= Ab((i+1):end,i)/Ab(i,i);
   %rep_Ab_row = repmat(Ab(i,:),n-i,1);
   %Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;
   Ab((i+1):end,i:end) = Ab((i+1):end,i:end)-l.*Ab(i,i:end);
end
A=Ab(:,1:n);
b=Ab(:,n+1);
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
   x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);
end