在 Python 中使用 Cramer 方法求解线性方程组答案

【问题标题】：Solving System of linear equation using Cramer's method in Python在 Python 中使用 Cramer 方法求解线性方程组
【发布时间】：2022-01-16 06:52:51
【问题描述】：

我正在尝试求解以下线性方程组：

10x1+ 40x2+ 70x3= 300
20x1+ 50x2+ 80x3= 360
30x1+ 60x2+ 80x3= 390

通过使用 Cramer 的方法从头实现一个函数：

def cramer(mat, constant): # takes the matrix and the costants
    
    D = np.linalg.det(mat) # calculating the determinant of the original matrix

# substitution of the column with costant and creating new matrix
    mat1 = np.array([constant, mat[1], mat[2]])
    mat2 = np.array([mat[0], constant, mat[2]])
    mat3 = np.array([mat[0], mat[1], constant])  
    
#calculatin determinant of the matrix
    D1 = np.linalg.det(mat1)
    D2 = np.linalg.det(mat2)
    D3 = np.linalg.det(mat3)
    
#finding the X1, X2, X3
    X1 = D1/D
    X2 = D2/D
    X3 = D3/D
    
    print(X1, X2, X3)

通过在系统上使用上述功能

a = np.array([[10, 40, 70],
             [20, 50, 80],
             [30, 60, 80]])

b = np.array([300, 360, 390])

我得到以下结果：

cramer(a,b)
-22.99999999999996 21.999999999999964 2.999999999999993

我已经使用 numpy 函数 np.linalg.solve 解决了系统问题，我得到了另一个结果：

x = np.linalg.solve(a, b)
[1. 2. 3.]

我无法在我写过的函数中发现公式错误。我应该在功能中进行哪些调整才能使其正常工作？

【问题讨论】：

标签： python numpy matrix linear-algebra linear-equation

【解决方案1】：

主要问题是如何选择a 的列，即您实际上是选择a 的行而不是列。您可以通过将矩阵初始化更改为此来修复它：

    mat1 = np.array([constant, mat[:,1], mat[:,2]])
    mat2 = np.array([mat[:,0], constant, mat[:,2]])
    mat3 = np.array([mat[:,0], mat[:,1], constant])

基本上mat[:,1] 说的是mat[all rows, column 1]。

【讨论】：

【解决方案2】：

TL;DR底部最优解。

要修复您当前的解决方案，您需要使用第二维并传递所有三个矩阵以一起计算行列式（这样您将获得稳定的浮点值）：

def cramer(mat, constant):
    D = np.linalg.det(mat)
    mat1 = np.array([constant, mat[:, 1], mat[:, 2]])
    mat2 = np.array([mat[:, 0], constant, mat[:, 2]])
    mat3 = np.array([mat[:, 0], mat[:, 1], constant])
    Dx = np.linalg.det([mat1, mat2, mat3])
    X = Dx/D
    print(X)

但是，您也不需要一一创建所有这些矩阵。相反，请使用下面描述的一系列numpy 操作。

首先，创建掩码，以便您可以使用它将a 中的值替换为b 中的值：

>>> mask = np.broadcast_to(np.diag([1,1,1]), [3, 3, 3]).swapaxes(0, 1)

array([[[1, 0, 0],
        [1, 0, 0],
        [1, 0, 0]],

       [[0, 1, 0],
        [0, 1, 0],
        [0, 1, 0]],

       [[0, 0, 1],
        [0, 0, 1],
        [0, 0, 1]]])

然后使用np.where得到三个矩阵，每个矩阵有一列替换为b：

>>> Ms = np.where(mask, np.repeat(b, 3).reshape(3, 3), a)

array([[[300,  40,  70],
        [360,  50,  80],
        [390,  60,  80]],

       [[ 10, 300,  70],
        [ 20, 360,  80],
        [ 30, 390,  80]],

       [[ 10,  40, 300],
        [ 20,  50, 360],
        [ 30,  60, 390]]])

然后，计算三个行列式并将a本身的行列式相除：

>>> np.linalg.det(Ms) / np.linalg.det(a)

array([1., 2., 3.])

把它们放在一起：

def cramer(a, b):
  mask = np.broadcast_to(np.diag([1,1,1]), [3, 3, 3]).swapaxes(0, 1)
  Ms = np.where(mask, np.repeat(b, 3).reshape(3, 3), a)
  return np.linalg.det(Ms) / np.linalg.det(a)

【讨论】：