在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

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)

Answer 2

主要问题是如何选择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]的东西。