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

Question

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

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.]

我无法在编写的函数中发现公式错误。功能上应该怎么调整才能正常使用？

Answer 1

主要问题是您如何 select a 的列，即您实际上是 selecting 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].

Answer 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 操作。

首先，创建掩码，这样您就可以使用它用 b:

中的值替换 a 中的值

>>> 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)

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

Solving System of linear equation using Cramer's method in Python

python

linear-equation

numpy

matrix

linear-algebra