Python 解线性方程组

线性方程组是各个方程的未知元的次数都是一次的方程组。解这样的方程组有两种方法：克拉默法则和矩阵消元法。

矩阵消元法

矩阵消元法。将线性方程组的增广矩阵通过行的初等变换化为行简化阶梯形矩阵 ，则以行简化阶梯形矩阵为增广矩阵的线性方程组与原方程组同解。当方程组有解时，将其中单位列向量对应的未知量取为非自由未知量，其余的未知量取为自由未知量，即可找出线性方程组的解。

这种方法适合手工解方程，通过编写程序来解方程这种方法基本行不通。因为行初等变换有 3 种：两行互换，某一行乘上一个非零实数，某一行乘上一个非零实数后加到另一行。看着很简单，但是关键问题是怎么让计算机知道哪两行互换？哪一行要乘上一个数？这个数是多少？

克拉默法则

因为上面的问题对于计算机来说是很难处理的，所以我们就换一种方法，这次使用克拉默法则。用克莱姆法则求解方程组有两个前提，一是方程的个数要等于未知量的个数，二是系数矩阵的行列式要不等于零。用克莱姆法则求解方程组实际上相当于用逆矩阵的方法求解线性方程组，它建立线性方程组的解与其系数和常数间的关系，但由于求解时要计算 n+1 个 n 阶行列式，其工作量常常很大，所以克莱姆法则常用于理论证明，很少用于具体求解。下面看一下克拉默法则的具体使用。

首先将方程组整理为如下形式：

a11*x1+a12*x2+...+a1n*xn=b1

a21*x1+a22*x2+...+a2n*xn=b2

...

an1*x1+an2*x2+...+ann*xn=bn

系数矩阵记为 A，将系数矩阵中的第 i 列换成对应的常数项，换好后的矩阵记为 Ai，那么 xi=|Ai|/|A|。下面我以 5 个未知数 5 个方程为例实现一下代码。

from numpy.linalg import det
from random import randint
from pprint import pprint
from copy import deepcopy
a = [[randint(1, 9)for i in range(5)]for j in range(5)]  # 系数矩阵
b = [randint(1, 9)for k in range(5)]  # 常数项
print('a=', end='')
pprint(a)
print(f'b={b}')
print()
# 克拉默法则
for i in range(5):
    a0 = deepcopy(a)  # 深拷贝，防止修改系数行列式
    for j in range(5):
        a0[j][i] = b[j]  # 将行列式中 xi 的系数变为对应的常数项
    print(f'x{i}={det(a0)/det(a)}')

其实还可以更简单，一个行列式等于提出某一行或某一列然后乘上对应的代数余子式做一个加和，在这里，我们将 |Ai| 里面的 bi 都提出来，得

|Ai|=∑bj*|Aji|（其中 |Aji| 为对应元素的代数余子式，1≤j≤n）将其代入

xi=|Ai|/|A| 得 xi=∑bj*|Aji|/|A|。由于表达式中有代数余子式，有行列式，为此我就可以想到 xi，A 的逆和 bi 有着某种联系，至于什么联系我们就先尝试把 A 的逆和 b 向量做矩阵乘法运算，得到一个向量，这个向量的第 i 个元素就是 xi 的值，既然如此，那么我就只要逆矩阵*常数向量就可以得出解向量 x 了，代码实现比上面那种方法简单太多了，一行代码就能求出解向量，代码如下：

# 系数矩阵的逆*常数向量
x = inv(a)@b
for i in range(5):
    print(f'x{i}={x[i]}')

最后附上完整的源代码和运行结果。

from numpy.linalg import det, inv
from random import randint
from pprint import pprint
from copy import deepcopy
a = [[randint(1, 9)for i in range(5)]for j in range(5)]  # 系数矩阵
b = [randint(1, 9)for k in range(5)]  # 常数项
print('a=', end='')
pprint(a)
print(f'b={b}')
print()
# 克拉默法则
for i in range(5):
    a0 = deepcopy(a)  # 深拷贝，防止修改系数行列式
    for j in range(5):
        a0[j][i] = b[j]  # 将行列式中 xi 的系数变为对应的常数项
    print(f'x{i}={det(a0)/det(a)}')
print()
# 系数矩阵的逆*常数向量
x = inv(a)@b
for i in range(5):
    print(f'x{i}={x[i]}')

可以发现，两种方法求出来的结果差不多，有一点误差是正常的，毕竟绝大多数浮点数计算机无法精确表示，要想精确表示小数可以使用模块 decimal。