解锁线性代数之门：用Python秒解任何线性方程组

解锁线性代数之门：用Python秒解任何线性方程组

在数学和计算机科学的世界中，解决线性方程组是一个基本而重要的任务。无论是在工程计算、数据分析还是科学研究中，线性方程组都扮演着不可或缺的角色。高斯消元法，以数学家卡尔·弗里德里希·高斯的名字命名，是解决这一问题的经典方法之一。它不仅在理论上具有重要意义，而且在实际应用中也极为高效。本文将带你深入探索高斯消元法的奥秘，从它的数学原理出发，一步步走向其Python实现。

要秒解任何线性方程组，我们首先要把一般形式的线性方程组转变成增广矩阵的形式，这样才方便用计算机程序运算。

一.把一般形式的线性方程组转变成增广矩阵

1.对于简单的线性方程组，最常见的方法是直接输入。比如说，我们看到题目如图，我们就看着这个图，在Python程序里面输入：

matrix=[[2,1,1,1],[6,2,1,-1],[-2,2,1,7]]

这样不就一下子把题目里的方程组变成标准的增广矩阵了吗？！

但是，如果我们要输入很多线性方程组，不仅自己看着眼花缭乱，还总是要去修改Python本体。

直接输入常量的程序维护性太差了！

2.我自己使用的方法是正则表达式。

正则表达式快速入门：https://www.zhengzeshi.com

我一开始这样应用正则表达式：

import re
# 示例字符串
s = "3a - 2b + 7c = 4"
# 正则表达式匹配位于任意英文字母前的数字，但结果不包含这些字母
pattern = r"-?\d*(?=[a-zA-Z])"
# 使用 re.findall 查找所有匹配项
matches = re.findall(pattern, s)
# 输出提取到的数字
print("提取到的数字包括:", numbers)
"""
提取到的数字包括: ['3', '2', '7']
"""

似乎这个程序是能够提取出来方程中的系数的。

但是，-?\d*(?=[a-zA-Z])这个正则表达式只能提取显式系数，而不适用于x+y+z=1这类包含隐式系数的方程，对于存在零系数的方程组，更是会忽略掉零系数！

那我们怎么解决隐式系数+1、-1和0的问题呢？

问一下AI怎么说~

这个Python程序先收集了方程组里一共有几个未知数变量，然后在方程组里面寻找各个变量，进行逻辑判断。

如果没有找到变量，则系数为0。

如果变量前面是空字符串，则系数为1。

如果变量前面是负号字符串，则系数为-1。

import numpy as np
import re

# 示例方程组
equations = ["x+y+z=1", "x+y=0"]

# 步骤1: 收集所有变量名
all_variables = set()
pattern = r"([a-zA-Z]+)"
for eq in equations:
    variables = re.findall(pattern, eq.split('=')[0])
    all_variables.update(variables)
all_variables = sorted(list(all_variables))  # 排序以保持一致性

# 步骤2: 构建方程系数矩阵
coeff_matrix = []
for eq in equations:
    coeff = [0] * len(all_variables)  # 初始化系数为0
    for var in all_variables:
        # 查找变量系数，如果没有找到，则系数为0
        match = re.search(rf"(-?\d*){var}", eq)
        if match:
            coeff_val = match.group(1)
            coeff[all_variables.index(var)] = int(coeff_val) if coeff_val not in ['', '-'] else -1 if coeff_val == '-' else 1
    coeff_matrix.append(coeff)

# 步骤3: 构建常数项向量
const_vector = [int(eq.split('=')[1]) for eq in equations]

# 输出结果
print("变量顺序:", all_variables)
print("系数矩阵:", coeff_matrix)
print("常数项向量:", const_vector)

"""
变量顺序: ['x', 'y', 'z']
系数矩阵: [[1, 1, 1], [1, 1, 0]]
常数项向量: [1, 0]
"""

3.对了，我还要给大家推荐一个很好用的专门的数学表达式解析库！——Sympy！

安装 SymPy: 如果您还没有安装 SymPy，可以通过 pip 来安装它：

pip install sympy

然后你就能用Sympy解析数学表达式并提取系数了:

from sympy import symbols, Eq, solve

# 定义符号
x = symbols('x')
# 定义方程 ax^2 + bx + c = 0
equation = Eq(2*x**2 + 3*x + 1, 0)
# 提取系数
a, b, c = equation.lhs.as_poly().coeffs()
print(f"系数a: {a}, 系数b: {b}, 系数c: {c}")

"""
系数a: 2, 系数b: 3, 系数c: 1
"""

分析a, b, c = equation.lhs.as_poly().coeffs()这句。

lhs 代表一个等式（Equation）对象的左侧部分（Left-Hand Side）。

让我们深入了解 SymPy 中的 as_poly 和 coeffs 方法，它们是处理多项式表达式时非常有用的工具。

as_poly 方法

as_poly 方法是 SymPy 表达式的一个方法，它将表达式转换为多项式对象。这个转换对于执行多项式特定的操作非常重要，比如提取系数、计算多项式的根，或者进行多项式的加减乘除。转换为多项式对象后，你可以利用 SymPy 提供的多项式特定功能，这些功能在一般的表达式上可能不可用。

基本用法:

expr.as_poly(*gens, **args)

• expr: 需要转换的 SymPy 表达式。
• *gens: 可选参数，用于指定表达式中作为多项式变量的符号。如果不指定，SymPy 会自动确定表达式中的变量。
• **args: 其他关键字参数，用于控制多项式的创建，例如指定域。

示例:

假设有表达式 $2x^2+3x+1$，转换为多项式对象：

from sympy import symbols
x = symbols('x')
expr = 2*x**2 + 3*x + 1
poly = expr.as_poly()

poly 现在是一个多项式对象，包含关于 x 的多项式 $2x^2+3x+1$。

coeffs 方法

coeffs 方法是多项式对象的一个方法，用于返回多项式的系数列表。这个方法在你需要分析多项式结构，例如提取系数来进行进一步计算或者验证时特别有用。

基本用法:

poly.coeffs()

这个方法返回一个列表，包含从最高次幂到最低次幂的系数。

示例:

继续上面的多项式 $2x^2+3x+1$的例子，使用 coeffs 方法提取系数：

coefficients = poly.coeffs()

coefficients 将是一个列表 [2, 3, 1]，对应于多项式 $2x^2+3x+1$ 的系数。

总结

• as_poly 方法允许你将 SymPy 表达式转换为多项式对象，这对于执行多项式特定的操作非常有用。
• coeffs 方法用于从多项式对象中提取系数列表，这可以帮助你分析多项式的结构或进行进一步的数学处理。

这两个方法在处理和分析多项式表达式时非常有用，特别是在数学、工程和科学计算中。

根据以上三个方法，相信你已经找到了系数矩阵和常量项。

下一步要做的就是根据矩阵来判断方程组解的情况。

要知道方程组有没有解，我们才能继续解方程啊！

二.计算矩阵秩RANK判断解的情况

1. 增广矩阵的秩等于系数矩阵的秩：这意味着所有的方程都是相容的，没有矛盾，因此至少存在一个解。

   • 如果这个共同的秩同时也等于变量的数量，那么系统有唯一解。
   • 如果这个共同的秩小于变量的数量，说明有不止一个解，系统有无穷多解，因为有自由变量存在。

2. 增广矩阵的秩大于系数矩阵的秩：这意味着方程组中至少有一个方程与其他方程不相容，导致系统无解。

   要怎么计算矩阵的秩呢？

   我们可以先把矩阵转换成行阶梯型，再计算非零行的数量。

我写了一个把普通矩阵转换到行阶梯形的函数ref：

#把增广矩阵转换成行阶梯形（REF）
def ref(A,b):
    n = len(A)
    # 前向消元
    for i in range(n):
        # 寻找主元
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k
        # 交换最大行至当前行
        A[i], A[maxRow] = A[maxRow], A[i]
        b[i], b[maxRow] = b[maxRow], b[i]
        # 将当前行的首元素归一化，并消去下方所有行的对应元素
        for k in range(i+1, n):
            c = -A[k][i] / A[i][i]
            for j in range(i, n):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]
            b[k] += c * b[i]
    return A, b

我们也可以不自己写方法来计算秩Rank，直接使用Numpy库的numpy.linalg.matrix_rank函数就能计算矩阵的秩Rank了。这个函数通过使用一个小的奇异值判断阈值来确定秩，它实际上是基于奇异值分解（SVD）来计算的。这种方法使得matrix_rank函数即便在处理数值上不稳定的矩阵时也能给出准确的秩。

如果想要更了解numpy.linalg.matrix_rank函数的具体实现，可以在以下网页的1973行到2088行找到numpy.linalg.matrix_rank的源代码。

https://github.com/numpy/numpy/blob/06b169bdcdc01d020c573d0af845327173e403c9/numpy/linalg/_linalg.py#L1973

A = coeff_matrix
b = const_vector
print('系数矩阵：',A)
print('常数项向量：',b)
#A, b = ref(A, b)#把增广矩阵转换成行阶梯形（REF）
augmented_matrix=np.column_stack((A,b))#合成增广矩阵
rank1 = np.linalg.matrix_rank(A)#求系数矩阵的秩
rank2=np.linalg.matrix_rank(augmented_matrix)#求增广矩阵的秩
n=len(all_variables)#未知数的个数
if rank1==rank2:
    if rank1==n:
        print('The system has a unique solution')#有唯一解
#调用高斯消元法求解方程组，输出结果
        x = gaussian_elimination(A, b)
        print('解为：',x)
    else:
        print('The system has infinite solutions')#有无穷多解
else:
    print('The system has no solution')#无解

三.高斯消元法解方程

高斯消元法的数学原理：https://cloud.tencent.com/developer/article/1087352

#高斯消元法求解线性方程组
def gaussian_elimination(A, b):
    n = len(A)
    # 前向消元
    for i in range(n):
        # 寻找主元
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k
        # 交换最大行至当前行
        A[i], A[maxRow] = A[maxRow], A[i]
        b[i], b[maxRow] = b[maxRow], b[i]
        # 将当前行的首元素归一化，并消去下方所有行的对应元素
        for k in range(i+1, n):
            c = -A[k][i] / A[i][i]
            for j in range(i, n):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]
            b[k] += c * b[i]
    # 回代
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = b[i] / A[i][i]
        for k in range(i-1, -1, -1):
            b[k] -= A[k][i] * x[i]
    return x

四.开箱即用的程序汇总

把之前谈及的方法全部组装起来，我们就能得到一个解任何线性方程组的Python程序了！

平时我们自己手动计算六元方程组都觉得很麻烦了，这个程序能解成千上万的方程组！

收藏加关注，直接在这里看全部源码。

import re
import numpy as np

#高斯消元法求解线性方程组
def gaussian_elimination(A, b):
    n = len(A)
    # 前向消元
    for i in range(n):
        # 寻找主元
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k
        # 交换最大行至当前行
        A[i], A[maxRow] = A[maxRow], A[i]
        b[i], b[maxRow] = b[maxRow], b[i]
        # 将当前行的首元素归一化，并消去下方所有行的对应元素
        for k in range(i+1, n):
            c = -A[k][i] / A[i][i]
            for j in range(i, n):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]
            b[k] += c * b[i]
    # 回代
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = b[i] / A[i][i]
        for k in range(i-1, -1, -1):
            b[k] -= A[k][i] * x[i]
    return x

#输入方程组
equations = []
T=int(input('How many equations do you have?\n'))
print('Enter the equation like 3x+4y+z=80:')
for i in range(T):
    equations.append(input())

# 步骤1: 收集所有变量名
all_variables = set()
pattern = r"([a-zA-Z]+)"
for eq in equations:
    variables = re.findall(pattern, eq.split('=')[0])
    all_variables.update(variables)
all_variables = sorted(list(all_variables)) # 排序以保持一致性

# 步骤2: 构建方程系数矩阵
coeff_matrix = []
for eq in equations:
    coeff = [0] * len(all_variables)  #初始化系数为0
    for var in all_variables:
        # 查找变量系数，如果没有找到，则系数为0
        match = re.search(rf"(-?\d*){var}", eq)
        if match:
            coeff_val = match.group(1)
            coeff[all_variables.index(var)] = int(coeff_val) if coeff_val not in ['', '-'] else -1 if coeff_val == '-' else 1
    coeff_matrix.append(coeff)

# 步骤3: 构建常数项向量
const_vector = [int(eq.split('=')[1]) for eq in equations]

# 步骤4: 根据行阶梯形REF的秩，判断方程组的解的情况
A = coeff_matrix
b = const_vector
print('系数矩阵：',A)
print('常数项向量：',b)
#A, b = ref(A, b)#把增广矩阵转换成行阶梯形（REF）
augmented_matrix=np.column_stack((A,b))#合成增广矩阵
rank1 = np.linalg.matrix_rank(A)#求系数矩阵的秩
rank2=np.linalg.matrix_rank(augmented_matrix)#求增广矩阵的秩
n=len(all_variables)#未知数的个数
if rank1==rank2:
    if rank1==n:
        print('The system has a unique solution')#有唯一解
#步骤5: 调用高斯消元法求解方程组，输出结果
        x = gaussian_elimination(A, b)
        print('解为：',x)
    else:
        print('The system has infinite solutions')#有无穷多解
else:
    print('The system has no solution')#无解