所谓的“线型系统”实际上可以理解为线性方程组,更严谨的定义为“包含未知数的项,未知数的次数只能为1的方程组”。研究线型系统实际上就是为了解出线性方程组的解,这一过程可以使用线型代数。
高斯消元法
将方程组中的未知数一个一个消去,这是最简单的思路,最终仅剩下一个未知数的时候进行求解,而后带入到其他方程中求得其余未知数的值。
下面总结了一般的代数中消元法求解方程组的三个基本操作,
- 一个方程左右两边同时乘以一个常数
- 一个方程加(减)另一个方程
- 交换两个方程的位置
通过上面的三种操作就可以改变线性系统达到消元的目的。
矩阵的方式可以很容易的表达一个线性系统,矩阵求解线型系统时候同样存在消元法,矩阵消元法求解将线性系统进一步简化,如下例,
{
x
+
2
y
+
4
z
=
7
3
x
+
7
y
+
2
z
=
−
11
2
x
+
3
y
+
3
z
=
1
\begin{cases}x+2y+4z=7\\3x+7y+2z=-11\\2x+3y+3z=1\end{cases}
⎩⎪⎨⎪⎧x+2y+4z=73x+7y+2z=−112x+3y+3z=1
表示为矩阵与向量相乘的形式,
(
1
2
4
3
7
2
2
3
3
)
⋅
(
x
y
z
)
=
(
7
−
11
1
)
\begin{pmatrix}1&2&4\\3&7&2\\2&3&3\end{pmatrix}\cdot \begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}7\\-11\\1\end{pmatrix}
⎝⎛132273423⎠⎞⋅⎝⎛xyz⎠⎞=⎝⎛7−111⎠⎞
转化为消元法形式,
(
1
2
4
7
3
7
2
−
11
2
3
3
1
)
\begin{pmatrix}1&2&4&7\\3&7&2&-11\\2&3&3&1\end{pmatrix}
⎝⎛1322734237−111⎠⎞
消元法中将未知数去掉是因为这些未知数并不能提供有用的信息,真正有价值的信息存在于系数矩阵中。
上面消元法形式的矩阵也称为“增广矩阵”。之前使用一般代数的消元法时的各步骤现在都可以在这个增广矩阵中进行。
{
x
+
2
y
+
4
z
=
7
3
x
+
7
y
+
2
z
=
−
11
2
x
+
3
y
+
3
z
=
1
⟹
{
x
+
2
y
+
4
z
=
7
y
−
10
z
=
−
32
−
y
−
5
z
=
−
13
⟹
{
x
+
2
y
+
4
z
=
7
y
−
10
z
=
−
32
z
=
3
\begin{cases}x+2y+4z=7\\3x+7y+2z=-11\\2x+3y+3z=1\end{cases}\Longrightarrow\begin{cases}x+2y+4z=7\\y-10z=-32\\-y-5z=-13\end{cases}\Longrightarrow\begin{cases}x+2y+4z=7\\y-10z=-32\\z=3\end{cases}
⎩⎪⎨⎪⎧x+2y+4z=73x+7y+2z=−112x+3y+3z=1⟹⎩⎪⎨⎪⎧x+2y+4z=7y−10z=−32−y−5z=−13⟹⎩⎪⎨⎪⎧x+2y+4z=7y−10z=−32z=3
与下方的增广矩阵一一对应,
(
1
2
4
7
3
7
2
−
11
2
3
3
1
)
⟹
(
1
2
4
7
0
1
−
10
−
32
0
−
1
−
5
−
13
)
⟹
(
1
2
4
7
0
1
−
10
−
32
0
0
−
15
−
45
)
\begin{pmatrix}1&2&4&7\\3&7&2&-11\\2&3&3&1\end{pmatrix}\Longrightarrow\begin{pmatrix}1&2&4&7\\0&1&-10&-32\\0&-1&-5&-13\end{pmatrix}\Longrightarrow\begin{pmatrix}1&2&4&7\\0&1&-10&-32\\0&0&-15&-45\end{pmatrix}
⎝⎛1322734237−111⎠⎞⟹⎝⎛10021−14−10−57−32−13⎠⎞⟹⎝⎛1002104−10−157−32−45⎠⎞
对应上面一般的代数问题中消元法对于方程组的操作,矩阵系统中的操作依次为,
- 矩阵的某一行乘以一个常数
- 矩阵的一行加(减)另一行
- 交换矩阵的两行
上面的消元过程可以得到如下的一个矩阵,
(
1
2
4
7
0
1
−
10
−
32
0
0
1
3
)
\begin{pmatrix}{\color{red} 1}&2&4&7\\0&{\color{red} 1}&-10&-32\\0&0&{\color{red} 1}&3\end{pmatrix}
⎝⎛1002104−1017−323⎠⎞
标红的三项称为系数矩阵的“主元”,消元的过程是使得下面一行比上面一行少一个未知数,相应的多出一个0。
对于一个n个未知数n个方程的线型系统而言,第 i i i行的主元实际上就是该行第 i i i个元素。上述消元的过程实际上可以视为从第一行开始,使用矩阵的三种操作,将主元的元素消减为1,这种方法也称为“高斯消元法”。
如果处理到某一行时,主元位置对应的元素为0,则选取该位置下面最大元素对应的一行与该行进行交换,这是为了减小计算误差(一种数值计算方法)。
高斯-约旦消元法(基础)
高斯消元法会得到主元位置均为1的增广矩阵,求解其余全部的未知数的解就需要使用“高斯-约旦消元法”(Gauss-Jordan elimination)。
前向过程(从上到下)
在这个方法中高斯消元法可以视为一个前向过程,还是以上一节中的线性系统为例,高斯消元法的前向过程得到的结果如下,
(
1
2
4
7
3
7
2
−
11
2
3
3
1
)
⟹
(
1
2
4
7
0
1
−
10
−
32
0
0
1
3
)
\begin{pmatrix}1&2&4&7\\3&7&2&-11\\2&3&3&1\end{pmatrix}\Longrightarrow\begin{pmatrix}{\color{red} 1}&2&4&7\\0&{\color{red} 1}&-10&-32\\0&0&{\color{red} 1}&3\end{pmatrix}
⎝⎛1322734237−111⎠⎞⟹⎝⎛1002104−1017−323⎠⎞
前向过程可以总结为,
- 选择最上面的主元(左上角),化为1。
- 主元下面所有行减去主元所在行的某个倍数,使得主元下面所有元素都为0。
- 最终第一行有n个非0元素,第二行n-1个非0元素,到最后一行只有一个非0元素,对应着最后一个未知数的值。
后向过程(从下到上)
约旦在此基础上添加了一个后向过程,与前向过程相反,后向过程是一个从下到上的过程,
- 选择最下面的主元
- 主元上面所有行减去主元所在行的某个倍数,使得主元上面所有元素都为0
由此得到整个线性系统中未知数的解。
代码实现
在实现Gauss-Jordan消元法之前首先对向量笔记中定义的Vector类添加一个方法,
def underlying_list(self):
"""返回向量对应的列表"""
return self._values.copy()
创建一个新的脚本 line_system.py,
# -*- coding:utf-8 -*-
from .matrix import Matrix
from .vector import Vector
class LinearSystem:
def __init__(self, A, b):
"""
:param A: 线性系统的矩阵,Matrix类对象
:param b: 线性系统最终的结果(每个方程的解),Vector类对象
"""
assert A.row_num() == len(b)
self._m = A.row_num()
self._n = A.col_num()
# 就目前而言,高斯-约旦法解决的是n个未知数对应n个方程的问题,以后回解决未知数与方程不等的问题
assert self._m == self._n # TODO: remove this restriction
# 创建增广矩阵
self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._m)]
def _max_row(self, index, n):
"""寻找主元下方最大元素所在的行并返回行号"""
best, ret = self.Ab[index][index], index
for i in range(index+1, n):
if self.Ab[i][index] > best:
best, ret = self.Ab[i][index], i
return ret
def _forward(self):
"""前向过程"""
n = self._m
# 首先遍历每行主元 self.Ab[i][i]如果为0,则与下方最大值的行交换
for i in range(n):
if self.Ab[i][i] == 0:
max_row = self._max_row(i, n)
self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]
# 主元归一
self.Ab[i] = self.Ab[i] /self.Ab[i][i]
# 主元一下所有的行减去主元的倍数进行消元
for j in range(i+1, n):
self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]
return None
def _backward(self):
"""后向过程"""
n = self._m
# self.Ab[i][i]为主元,将其上方所有元素变为0,注意此时从最后一行开始
for i in range(n-1, -1, -1):
for j in range(i-1, -1, -1):
self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]
return None
def gauss_jordan_elimination(self):
"""高斯-约旦消元法"""
self._forward()
self._backward()
def fancy_print(self):
for i in range(self._m):
print(" ".join(str(self.Ab[i][j]) for j in range(self._n)), end=" ")
print("|", self.Ab[i][-1])
return None
对上述代码进行测试,
from playLA.matrix import Matrix
from playLA.vector import Vector
from playLA.linear_system import LiearSystem
if __name__ == "__main__":
A = Matrix([[1, 2, 4], [3, 7, 2], [2, 3, 3]])
b = Vector([7, -11, 1])
ls = LinearSystem(A, b)
ls.gauss_jordan_elimination()
ls.fancy_print()
# 打印结果
# 1.0 0.0 0.0 | -1.0
# 0.0 1.0 0.0 | -2.0
# 0.0 0.0 1.0 | 3.0
目前的局限
目前的消元法还是有局限,这是之后需要解决的,具体可以总结如下,
- 目前仅支持n个未知数对应n个方程求解
- 线性矩阵只能是唯一解的
而实际的线性系统中不一定未知数数量等于方程数量。其次线性系统可能是无解的,或者是具有无穷解的。