I am trying to solve the following system of linear equations:

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

by using Cramer's method implementing a function by scratch:

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)

By using the above function on the system

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

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

I get the following result:

cramer(a,b)
-22.99999999999996 21.999999999999964 2.999999999999993

I have solved the system using the numpy function np.linalg.solve and I get another result:

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

I cannot spot the formula error in the function I have witten. What should I adjust in the fuction in order to make it working properly?

CodePudding user response：

The main problem is how you select the columns of a, i.e. you are actually selecting the rows of a rather than the columns. You can fix it by changing the matrix initializations to this:

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

Basically mat[:,1] is saying something like mat[all rows, column 1].

CodePudding user response：

TL;DR Optimal solution at the bottom.

To fix your current solution you need to use the second dimensions and pass all three matrices to compute determinants together (this way you will get stable floating point values):

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)

However, you don't need to create all these matrices one by one either. Instead, use a series of numpy manipulations described below.

First, create the mask to so you can then use it to replace values in a by values from 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]]])

Then use np.where to get three matrices, each with one column replaced by 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]]])

Then, compute three determinants and divide the determinant of a itself:

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

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

Putting it all together:

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)