Search code examples
pythonarraysnumpymultiplication

Numpy - Find product Ax=b of custom A nxn matrix and B nx1


Want to find x , from Ax=b .First of all I have declared two matrices , A which is nxn and B nx1 . Their formulas can be seen below. For A: enter image description here

and for b : enter image description here

The matrices can take any n. In my code I'm giving it a value of 10. I declare them by set them to zero first . Then f I declare each element of 1,2 and n-1 and n line for both matrices . For A I also loop each number to get the desired look , each same number goes one col and row forward , from n=2 until n-2. Then for calculating the Ax=b , for finding x , I use the multiplication :

x = np.dot(np.linalg.inv(A), b) 

but doesn't get the correct answer . Any help ? I'm posting also my code in its entirety:

import numpy as np

n = 10

################## AAAAA matrix #############################################
A = np.zeros([n, n], dtype=float)  # initialize to f zeros

# ------------------first row
A[0][0] = 6
A[0][1] = -4
A[0][2] = 1
# ------------------second row
A[1][0] = -4
A[1][1] = 6
A[1][2] = -4
A[1][3] = 1
# --------------two last rows-----
# n-2 row
A[- 2][- 1] = -4
A[- 2][- 2] = 6
A[- 2][- 3] = -4
A[- 2][- 4] = 1
# n-1 row
A[- 1][- 1] = 6
A[- 1][- 2] = -4
A[- 1][- 3] = 1

# --------------------------- from second to n-2 row --------------------------#
j = 0
for i in range(2, n - 2):
    if j == (n - 4):
        break
    A[i][j] = 1
    j = j + 1

j = 1
for i in range(2, n - 2):
    if j == (n - 3):
        break
    A[i][j] = -4
    j = j + 1

j = 2
for i in range(2, n - 2):
    if j == (n - 2):
        break
    A[i][j] = 6
    j = j + 1

j = 3
for i in range(2, n - 2):
    if j == (n - 1):
        break
    A[i][j] = -4
    j = j + 1

j = 4
for i in range(2, n - 2):
    if j == (n):
        break
    A[i][j] = 1
    j = j + 1
# -----------------------------end coding of 2nd to n-2 r-------------#
print("\nMatrix A is : \n", A)

####### b matrix ######################################
b = np.zeros(n,float).reshape((n,1))
b[0] = 3
b[1] = -1
#b[len(b) - 1] = 3
#b[len(b) - 2] = -1
b[[0,-1]]=3; b[[1,-2]]=-1

print("\nMatrix b is \n", b)

#################### result ########################
x = np.dot(np.linalg.inv(A), b)

print("\n\n The result is : \n", x)

Actually I get all 1s for result , as you can see :

enter image description here


Solution

  • Using matmul method makes it:

    import numpy as np
    A = np.array([[1,2], [3,4]], dtype=float)
    B = np.array([[3], [4]])
    print(np.matmul(np.linalg.inv(A), B))