带有 Sympy 和 Scipy/Numpy 的矩阵的三角系统

Question

我刚刚开始自学线性代数，并且正在尝试自学一些基本的 Python 模块。如果可以的话，我决定学习 Scipy 和 Numpy 的线性代数组件，以及一些 Sympy。

我遇到过以下练习：

我已经手动解决了第一步（我相信）：

我想使用 Numpy 和/或 Scipy 以及 Sympy 来解决这个问题。我今天花了一些时间试图弄清楚使用什么结构来解决这个问题（显然是 Numpy 数组，然后是 Sympy 矩阵），但我不确定如何将其转换为三角形数组。我在求解方面见过很多，但没有什么能简单地产生上三角系统。

是否有人对使用其中一个或所有这些包执行此操作有任何建议？

Answer 1

除上述内容外，我还注意到here 关于对矩阵进行基本运算。这样做似乎总是很尴尬，你可能会写一些东西来帮助你在学习如何使用矩阵时更自然地工作。这是一个do 例程，它可能是一个帮助或如何编写自己的建议：

def do(self, s, i, j, k=1):
    """return self after doing row (`s` starts with 'r') or column (`s` starts with 'c') operation.

    Examples
    ========

    >>> from sympy import Matrix
    >>> Matrix.do = do
    >>> m = Matrix(3, 3, list(range(9))); m
    Matrix([
    [0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]])

    swap rows 0 and 1

    >>> m.do('rs',0,1)
    Matrix([
    [3, 4, 5],
    [0, 1, 2],
    [6, 7, 8]])

    To row 0, add row 1 multiplied by 2

    >>> m.do('ra',0,1,2)
    Matrix([
    [6, 9, 12],
    [3, 4,  5],
    [6, 7,  8]])

    Multiply row 0 by 2

    >>> m.do('rm',0,2)
    Matrix([
    [0, 2, 4],
    [3, 4, 5],
    [6, 7, 8]])

    To operate on columns, replace the 'r' with 'c' in the
    examples above.

    If you want to keep working with the modified matrix
    you can either re-assign the output to the variable used for
    the original matrix (e.g. `m = m.do('rs', 0, 1)`) or make the
    modifications in place by appending an 'i' to `s`:

    >>> m.do('rsi',0,1)
    >>> m
    Matrix([
    [3, 4, 5],
    [0, 1, 2],
    [6, 7, 8]])
    >>> m.do('csi',0,1)
    >>> m
    Matrix([
    [4, 3, 5],
    [1, 0, 2],
    [7, 6, 8]])
    """
    rc, op = s.lower()[:2]
    assert op in 'sma'
    assert rc in 'rc'
    if len(s) > 2:
        # in place
        if rc == 'r':
            if op == 's':
                self.row_swap(i,j)
            elif op == 'm':
                self.row_op(i, lambda v, _: v*j)
            else:
                self.row_op(i, lambda v, _: v + self[j, _]*k)
        else:
            if op == 's':
                self.col_swap(i,j)
            elif op == 'm':
                self.col_op(i, lambda v, _: v*j)
            else:
                self.col_op(i, lambda v, _: v + self[_, j]*k)
    else:
        f = getattr(self, 'elementary_%s_op' % (
            'row' if rc=='r' else 'col'))
        if op == 's':
            return f('n<->m', None, None, i, j)
        elif op == 'a':
            return f('n->n+km', None, k, i, j)
        else:
            return f('n->kn', i, j)

Answer 2

In [731]: arr = np.array([[2,-3,0,3],[4,-5,1,7],[2,-1,-3,5]])

您的两个步骤可以通过以下方式完成：

减去第一行的倍数：

In [732]: arr -= arr[0,:]*np.array([0,2,1])[:,None]
In [733]: arr
Out[733]: 
array([[ 2, -3,  0,  3],
       [ 0,  1,  1,  1],
       [ 0,  2, -3,  2]])

减去第二行的倍数：

In [734]: arr -= arr[1,:]*np.array([0,0,2])[:,None]
In [735]: arr
Out[735]: 
array([[ 2, -3,  0,  3],
       [ 0,  1,  1,  1],
       [ 0,  0, -5,  0]])

但我不确定这是否有助于您学习线性代数。虽然我想出了如何做到这一点，但这并不是我们在使用 numpy 时通常会做的事情。

Answer 3

我不确定您到底想做什么，但这是您可以在 SymPy 中执行一些基本线性代数例程的方法：

In [1]: import sympy as sym

In [2]: M = sym.Matrix([[2,-3,0,3],[4,-5,1,7],[2,-1,-3,5]])

In [3]: M
Out[3]: 
⎡2  -3  0   3⎤
⎢            ⎥
⎢4  -5  1   7⎥
⎢            ⎥
⎣2  -1  -3  5⎦

In [4]: M.echelon_form()
Out[4]: 
⎡2  -3   0   3⎤
⎢             ⎥
⎢0  2    2   2⎥
⎢             ⎥
⎣0  0   -20  0⎦

In [5]: M.rref()
Out[5]: 
⎛⎡1  0  0  3⎤           ⎞
⎜⎢          ⎥           ⎟
⎜⎢0  1  0  1⎥, (0, 1, 2)⎟
⎜⎢          ⎥           ⎟
⎝⎣0  0  1  0⎦           ⎠

In [9]: sym.linsolve(M)
Out[9]: {(3, 1, 0)}