高斯消元法
设线性方程组
A
x
=
b
Ax=b
Ax=b，其中
(
A
(
1
)
,
b
(
1
)
)
=
(
a
11
(
1
)
a
12
(
1
)
a
13
(
1
)
…
a
1
n
(
1
)
b
1
(
1
)
a
21
(
1
)
a
22
(
1
)
a
23
(
1
)
…
a
2
n
(
1
)
b
2
(
1
)
a
31
(
1
)
a
32
(
1
)
a
33
(
1
)
…
a
3
n
(
1
)
b
3
(
1
)
⋮
⋮
⋮
⋱
⋮
⋮
a
n
1
(
1
)
a
n
2
(
1
)
a
n
3
(
1
)
…
a
n
n
(
1
)
b
n
(
1
)
)
(A^{(1)},b^{(1)})= \left(\begin{array}{cccccc} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ a_{21}^{(1)} & a_{22}^{(1)} & a_{23}^{(1)} & \ldots & a_{2 n}^{(1)} & b_{2}^{(1)} \\ a_{31}^{(1)} & a_{32}^{(1)} & a_{33}^{(1)} & \ldots & a_{3 n}^{(1)} & b_{3}^{(1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n 1}^{(1)} & a_{n 2}^{(1)} & a_{n 3}^{(1)} & \ldots & a_{n n}^{(1)} & b_{n}^{(1)} \end{array}\right)
(A(1),b(1))=⎝⎜⎜⎜⎜⎜⎜⎛a11(1)a21(1)a31(1)⋮an1(1)a12(1)a22(1)a32(1)⋮an2(1)a13(1)a23(1)a33(1)⋮an3(1)………⋱…a1n(1)a2n(1)a3n(1)⋮ann(1)b1(1)b2(1)b3(1)⋮bn(1)⎠⎟⎟⎟⎟⎟⎟⎞
用高斯消元法时，步骤如下：
- 计算
−
l
i
1
=
−
a
i
1
(
1
)
a
11
(
1
)
,
(
i
=
2
,
3
,
.
.
.
,
n
)
-l_{i1}=-\frac{a^{(1)}_{i1}}{a^{(1)}_{11}},(i=2,3,...,n)
−li1=−a11(1)ai1(1),(i=2,3,...,n)，然后将矩阵的第一行乘这个系数分别加到每一行，这样，第一列的对角元一下就全是0。
( A ( 2 ) , b ( 2 ) ) = ( a 11 ( 1 ) a 12 ( 1 ) a 13 ( 1 ) … a 1 n ( 1 ) b 1 ( 1 ) 0 a 22 ( 2 ) a 23 ( 2 ) … a 2 n ( 2 ) b 2 ( 2 ) 0 a 32 ( 2 ) a 33 ( 2 ) … a 3 n ( 2 ) b 3 ( 2 ) ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 a n 2 ( 2 ) a n 3 ( 2 ) … a n n ( 2 ) b n ( 2 ) ] \left(A^{(2)}, b^{(2)}\right)=\left(\begin{array}{cccccc} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ 0 & a_{22}^{(2)} & a_{23}^{(2)} & \ldots & a_{2 n}^{(2)} & b_{2}^{(2)} \\ 0 & a_{32}^{(2)} & a_{33}^{(2)} & \ldots & a_{3 n}^{(2)} & b_{3}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & a_{n 2}^{(2)} & a_{n 3}^{(2)} & \ldots & a_{n n}^{(2)} & b_{n}^{(2)} \end{array}\right] (A(2),b(2))=⎝⎜⎜⎜⎜⎜⎜⎛a11(1)00⋮0a12(1)a22(2)a32(2)⋮an2(2)a13(1)a23(2)a33(2)⋮an3(2)………⋱…a1n(1)a2n(2)a3n(2)⋮ann(2)b1(1)b2(2)b3(2)⋮bn(2)⎦⎥⎥⎥⎥⎥⎥⎤
然后用除去第一行和第一列的矩阵重复上述步骤，n-1步之后，将会得到
( A ( n ) , b ( n ) ) = ( a 11 ( 1 ) a 12 ( 1 ) a 13 ( 1 ) … a 1 n ( 1 ) b 1 ( 1 ) a 22 ( 2 ) a 23 ( 2 ) … a 2 n ( 2 ) b 2 ( 2 ) a 33 ( 3 ) … a 3 n ( 3 ) b 3 ( 3 ) ⋱ ⋮ ⋮ a n n ( n ) b n ( n ) ) \left(A^{(n)}, b^{(n)}\right)=\left(\begin{array}{llllll} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ & a_{22}^{(2)} & a_{23}^{(2)} & \ldots & a_{2 n}^{(2)} & b_{2}^{(2)} \\ & & a_{33}^{(3)} & \ldots & a_{3 n}^{(3)} & b_{3}^{(3)} \\ & & & \ddots & \vdots & \vdots \\ & & & & a_{n n}^{(n)} & b_{n}^{(n)} \end{array}\right) (A(n),b(n))=⎝⎜⎜⎜⎜⎜⎜⎛a11(1)a12(1)a22(2)a13(1)a23(2)a33(3)………⋱a1n(1)a2n(2)a3n(3)⋮ann(n)b1(1)b2(2)b3(3)⋮bn(n)⎠⎟⎟⎟⎟⎟⎟⎞
这个时候消元完成，然后进行求解，如下：
{
x
n
=
b
n
(
n
)
/
a
n
n
(
n
)
x
i
=
(
b
i
(
i
)
−
∑
j
=
i
+
1
n
a
i
j
(
i
)
x
j
)
/
a
i
i
(
i
)
,
i
=
n
−
1
,
n
−
2
,
⋯
,
1
\left\{\begin{array}{l} x_{n}=b_{n}^{(n)} / a_{n n}^{(n)} \\ x_{i}=\left(b_{i}^{(i)}-\sum_{j=i+1}^{n} a_{i j}^{(i)} x_{j}\right) / a_{i i}^{(i)}, \quad i=n-1, n-2, \cdots, 1 \end{array}\right.
{xn=bn(n)/ann(n)xi=(bi(i)−∑j=i+1naij(i)xj)/aii(i),i=n−1,n−2,⋯,1
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
)
≠
0
a^{(i)}_{i,i}\neq0
ai,i(i)=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