代数解矩阵方程¶

使用 SymPy 解矩阵（线性）方程。例如，解 \( \left[\begin{array}{cc} c & d\\1 & -e\end{array}\right] \left[\begin{array}{cc} x\\y\end{array}\right] = \left[\begin{array}{cc} 2\\0\end{array}\right] \) 得到 \( \left[\begin{array}{cc} x\\y\end{array}\right] = \left[\begin{array}{cc} \frac{2e}{ce+d}\\\frac{2}{ce+d}\end{array}\right]\).

要考虑的替代方案¶

如果您的矩阵和常量向量仅包含数字，而不是符号，例如 \(\left[\begin{array}{cc} 1 & 2\\3 & 4\end{array}\right] \left[\begin{array}{cc} x\\y\end{array}\right] = \left[\begin{array}{cc} 2\\0\end{array}\right]\)，您可以使用以下其他免费开源软件包来代替 SymPy
- NumPy 的 numpy.linalg.solve()
- SciPy 的 scipy.linalg.solve()
- mpmath 的 lu_solve()
求解矩阵方程等效于求解线性方程组，因此如果您愿意，您可以代数解方程组
如果你将你的问题表述为线性方程组，并希望将其转换为矩阵形式，可以使用 linear_eq_to_matrix()，然后按照本指南中的步骤操作。

求解矩阵方程¶

以下是用 SymPy 的 sympy.matrices.matrixbase.MatrixBase.solve() 求解矩阵方程的示例。我们使用标准的矩阵方程公式 \(Ax=b\)，其中

\(A\) 是表示线性方程中系数的矩阵
\(x\) 是要解的未知数列向量
\(b\) 是常数列向量，其中每一行都是一个方程的值

>>> from sympy import init_printing
>>> init_printing(use_unicode=True)

>>> from sympy import symbols
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> A
⎡c  d ⎤
⎢     ⎥
⎣1  -e⎦
>>> b = Matrix([2, 0])
>>> b
⎡2⎤
⎢ ⎥
⎣0⎦
>>> A.solve(b)
⎡  2⋅e  ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢       ⎥
⎢   2   ⎥
⎢───────⎥
⎣c⋅e + d⎦

指南¶

矩阵通常必须是方阵¶

矩阵 \(A\) 通常必须是方阵，才能表示一个未知数数量与方程数量相同的线性方程组。如果不是，SymPy 将给出错误 ShapeError: `self` and `rhs` must have the same number of rows.

矩阵必须是方阵的要求的例外来自 SymPy 对 Moore-Penrose pseudoinverse 的使用。

求解矩阵方程的方法¶

SymPy 的矩阵求解方法 sympy.matrices.matrixbase.MatrixBase.solve() 可以使用几种不同的方法，这些方法在 API 参考链接中列出。根据矩阵的性质，特定方法可能更有效。默认情况下，将使用高斯-约旦消元法。

在 solve 中指定方法等同于使用专门的求解函数。例如，使用 solve 且 method='LU' 等同于调用 LUsolve()。

用同一个矩阵求解多个矩阵方程¶

如果你需要用同一个矩阵 \(A\) 但不同的常数向量 \(b\) 重复求解矩阵方程，则使用以下方法之一效率更高。

你可以通过 LUsolve() 使用 LU 分解

>>> from sympy import symbols, Matrix, eye, simplify
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> A
⎡c  d ⎤
⎢     ⎥
⎣1  -e⎦
>>> b = Matrix([2, 0])
>>> b
    ⎡2⎤
    ⎢ ⎥
    ⎣0⎦
>>> solution = A.LUsolve(b)
>>> solution
    ⎡  2⋅e  ⎤
    ⎢───────⎥
    ⎢c⋅e + d⎥
    ⎢       ⎥
    ⎢   2   ⎥
    ⎢───────⎥
    ⎣c⋅e + d⎦
>>> # Demonstrate that solution is correct
>>> simplify(A * solution)
    ⎡2⎤
    ⎢ ⎥
    ⎣0⎦
>>> b2 = Matrix([4, 0])
>>> b2
    ⎡4⎤
    ⎢ ⎥
    ⎣0⎦
>>> solution2 = A.LUsolve(b2)
>>> solution2
    ⎡  4⋅e  ⎤
    ⎢───────⎥
    ⎢c⋅e + d⎥
    ⎢       ⎥
    ⎢   4   ⎥
    ⎢───────⎥
    ⎣c⋅e + d⎦
>>> # Demonstrate that solution2 is correct
>>> simplify(A * solution2)
    ⎡4⎤
    ⎢ ⎥
    ⎣0⎦

另一种方法是计算逆矩阵，但这几乎总是更慢，对于更大的矩阵来说要慢得多。如果计算效率不是优先考虑，你可以使用 inv()

>>> from sympy import symbols, Matrix, simplify
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> b = Matrix([2, 0])
>>> b
    ⎡2⎤
    ⎢ ⎥
    ⎣0⎦
>>> b2 = Matrix([4, 0])
>>> b2
    ⎡4⎤
    ⎢ ⎥
    ⎣0⎦
>>> inv = A.inv()
>>> inv
    ⎡   e        d   ⎤
    ⎢───────  ───────⎥
    ⎢c⋅e + d  c⋅e + d⎥
    ⎢                ⎥
    ⎢   1       -c   ⎥
    ⎢───────  ───────⎥
    ⎣c⋅e + d  c⋅e + d⎦
>>> # Solves Ax = b for x
>>> solution = inv * b
>>> solution
    ⎡  2⋅e  ⎤
    ⎢───────⎥
    ⎢c⋅e + d⎥
    ⎢       ⎥
    ⎢   2   ⎥
    ⎢───────⎥
    ⎣c⋅e + d⎦
>>> # Demonstrate that solution is correct
>>> simplify(A * solution)
    ⎡2⎤
    ⎢ ⎥
    ⎣0⎦
>>> # Solves Ax = b2 for x
>>> solution2 = inv * b2
>>> solution2
    ⎡  4⋅e  ⎤
    ⎢───────⎥
    ⎢c⋅e + d⎥
    ⎢       ⎥
    ⎢   4   ⎥
    ⎢───────⎥
    ⎣c⋅e + d⎦
>>> # Demonstrate that solution2 is correct
>>> simplify(A * solution2)
    ⎡4⎤
    ⎢ ⎥
    ⎣0⎦

确定大型符号矩阵的逆矩阵在计算上可能不可行。

使用符号矩阵¶

操纵符号矩阵的计算复杂度会随着矩阵大小的增加而迅速增加。例如，符号矩阵行列式的项数随着矩阵维度的阶乘而增加。因此，可以求解的矩阵的最大维数比数值矩阵更有限。例如，这个 4x4 符号矩阵的行列式有 24 项，每项有四个元素

>>> from sympy import MatrixSymbol
>>> A = MatrixSymbol('A', 4, 4).as_explicit()
>>> A
⎡A₀₀  A₀₁  A₀₂  A₀₃⎤
⎢                  ⎥
⎢A₁₀  A₁₁  A₁₂  A₁₃⎥
⎢                  ⎥
⎢A₂₀  A₂₁  A₂₂  A₂₃⎥
⎢                  ⎥
⎣A₃₀  A₃₁  A₃₂  A₃₃⎦
>>> A.det()
A₀₀⋅A₁₁⋅A₂₂⋅A₃₃ - A₀₀⋅A₁₁⋅A₂₃⋅A₃₂ - A₀₀⋅A₁₂⋅A₂₁⋅A₃₃ + A₀₀⋅A₁₂⋅A₂₃⋅A₃₁ +
A₀₀⋅A₁₃⋅A₂₁⋅A₃₂ - A₀₀⋅A₁₃⋅A₂₂⋅A₃₁ - A₀₁⋅A₁₀⋅A₂₂⋅A₃₃ + A₀₁⋅A₁₀⋅A₂₃⋅A₃₂ +
A₀₁⋅A₁₂⋅A₂₀⋅A₃₃ - A₀₁⋅A₁₂⋅A₂₃⋅A₃₀ - A₀₁⋅A₁₃⋅A₂₀⋅A₃₂ + A₀₁⋅A₁₃⋅A₂₂⋅A₃₀ +
A₀₂⋅A₁₀⋅A₂₁⋅A₃₃ - A₀₂⋅A₁₀⋅A₂₃⋅A₃₁ - A₀₂⋅A₁₁⋅A₂₀⋅A₃₃ + A₀₂⋅A₁₁⋅A₂₃⋅A₃₀ +
A₀₂⋅A₁₃⋅A₂₀⋅A₃₁ - A₀₂⋅A₁₃⋅A₂₁⋅A₃₀ - A₀₃⋅A₁₀⋅A₂₁⋅A₃₂ + A₀₃⋅A₁₀⋅A₂₂⋅A₃₁ +
A₀₃⋅A₁₁⋅A₂₀⋅A₃₂ - A₀₃⋅A₁₁⋅A₂₂⋅A₃₀ - A₀₃⋅A₁₂⋅A₂₀⋅A₃₁ + A₀₃⋅A₁₂⋅A₂₁⋅A₃₀

求解它的矩阵方程大约需要一分钟，而类似的 3x3 矩阵则不到一秒钟。矩阵中无关的符号条目越多，操纵它就越有可能变慢。这个例子，即找到所有元素都是独立符号的矩阵的通解，是极端情况，因此对于其大小的矩阵来说是最慢的。

加速求解矩阵方程¶

以下是一些建议

如果矩阵元素为零，请确保它们被识别为零。你可以通过将它们设为零或应用假设来做到这一点。
选择适合矩阵属性的求解方法，例如厄米特、对称或三角形。参考求解矩阵方程的方法。
使用 DomainMatrix 类，它操作起来可能更快，因为它限制了矩阵元素的域。

使用求解结果¶

将求解结果用作向量¶

你可以将求解结果用作向量。例如，要证明求解结果 \(x\) 是正确的，可以将其乘以矩阵 \(A\) 并验证它是否产生了常数向量 \(b\)

>>> from sympy import symbols, simplify
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> b = Matrix([2, 0])
>>> solution = A.solve(b)
>>> solution
    ⎡  2⋅e  ⎤
    ⎢───────⎥
    ⎢c⋅e + d⎥
    ⎢       ⎥
    ⎢   2   ⎥
    ⎢───────⎥
    ⎣c⋅e + d⎦
>>> # Not immediately obvious whether this result is a zeroes vector
>>> (A * solution) - b
    ⎡ 2⋅c⋅e      2⋅d      ⎤
    ⎢─────── + ─────── - 2⎥
    ⎢c⋅e + d   c⋅e + d    ⎥
    ⎢                     ⎥
    ⎣          0          ⎦
>>> # simplify reveals that this result is a zeroes vector
>>> simplify((A * solution) - b)
    ⎡0⎤
    ⎢ ⎥
    ⎣0⎦

请注意，我们必须使用 simplify() 来使 SymPy 简化矩阵元素中的表达式，使其立即显而易见，求解结果是正确的。

从求解结果中提取元素¶

因为你可以在列向量中遍历元素，所以你可以使用标准的 Python 技术来提取它的元素。例如，你可以使用列表推导创建一个包含元素的列表

>>> [element for element in solution]
    ⎡  2⋅e       2   ⎤
    ⎢───────, ───────⎥
    ⎣c⋅e + d  c⋅e + d⎦

或者你可以通过下标提取单个元素

>>> solution[0]
      2⋅e
    ───────
    c⋅e + d

无解的方程¶

如果矩阵的行列式为零，则用它表示的矩阵方程无解

>>> from sympy import symbols
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c*e**2, d*e], [c*e, d]])
>>> A
    ⎡   2     ⎤
    ⎢c⋅e   d⋅e⎥
    ⎢         ⎥
    ⎣c⋅e    d ⎦
>>> b = Matrix([2, 0])
>>> A.LUsolve(b)
Traceback (most recent call last):
    ...
NonInvertibleMatrixError: Matrix det == 0; not invertible.

报告错误¶

如果你在矩阵求解函数中发现错误，请在 SymPy 邮件列表上发布问题。在问题解决之前，你可以使用其他方法中列出的其他方法。