Could not find a correct inverse value of a matrix by using numpy

Question

I tried to use numpy to solve a matrix inverse, but could not get it to work. The Python code is shown below:

import copy
import numpy as np
from numpy.linalg import *

K=[[0.13535533905932737, -0.03535533905932737, 0.0, 0.0, -0.1],       
   [-0.03535533905932737, 0.13535533905932737, 0.0, -0.1, 0.0], 
   [0.0, 0.0, 0.13535533905932737, 0.03535533905932737, -0.03535533905932737], 
   [0.0, -0.1, 0.03535533905932737, 0.13535533905932737, -0.03535533905932737],
   [-0.1, 0.0, -0.03535533905932737, -0.03535533905932737, 0.13535533905932737]]

newK = np.array(K)
newK = inv(newK)

And the result is shown as:

[[ -2.74615429e+16  -2.74615429e+16   1.89368970e+00  -2.74615429e+16
   -2.74615429e+16]
 [ -2.74615429e+16  -2.74615429e+16  -1.37357926e+00  -2.74615429e+16
   -2.74615429e+16]
 [  9.23495156e-01  -1.84699031e+00   8.84484598e+00  -3.69398063e+00
    1.84699031e+00]
 [ -2.74615429e+16  -2.74615429e+16  -2.52873329e+00  -2.74615429e+16
   -2.74615429e+16]
 [ -2.74615429e+16  -2.74615429e+16   3.04884372e+00  -2.74615429e+16
   -2.74615429e+16]]

I also tried using MATLAB to check this value, but it turns out that:

    a =

    0.1340   -0.0350         0         0   -0.1000
   -0.0350    0.1350         0   -0.1000         0
         0         0    0.1350    0.0350   -0.0350
         0   -0.1000    0.0350    0.1350   -0.0350
   -0.1000         0   -0.0350   -0.0350    0.1350

>> inv(a)

ans =

   1.0e+03 *

   -1.0000   -1.0000    0.0000   -1.0000   -1.0000
   -1.0000   -0.9808   -0.0033   -0.9841   -0.9967
    0.0000   -0.0033    0.0089   -0.0044    0.0011
   -1.0000   -0.9841   -0.0044   -0.9785   -0.9956
   -1.0000   -0.9967    0.0011   -0.9956   -0.9911

Could someone please help me with this?

Answer 1

Firstly, the two matrices for the two programs are not identical:

numpy:

In [295]: np.linalg.det(newK)
Out[295]: -4.5051744884111939e-21

matlab:

>> det(a)
ans =
    -1.067499999999957e-07

Using the same matrix in numpy:

In [296]: K2=np.array([[0.1340, -0.0350, 0, 0, -0.1000],[-0.0350, 0.1350, 0, -0.1000, 0],[0, 0, 0.1350, 0.0350, -0.0350],[0, -0.1000, 0.0350, 0.1350, -0.0350],[-0.1000, 0, -0.0350, -0.0350, 0.135]])

In [297]: np.linalg.det(K2)
Out[297]: -1.0674999999999714e-07

And the inverse comes out right:

In [298]: np.linalg.inv(K2)
Out[298]: 
array([[ -1.00000000e+03,  -1.00000000e+03,  -2.93090593e-14,
         -1.00000000e+03,  -1.00000000e+03],
       [ -1.00000000e+03,  -9.80796253e+02,  -3.27868852e+00,
         -9.84074941e+02,  -9.96721311e+02],
       [ -5.26327952e-14,  -3.27868852e+00,   8.85245902e+00,
         -4.42622951e+00,   1.14754098e+00],
       [ -1.00000000e+03,  -9.84074941e+02,  -4.42622951e+00,
         -9.78501171e+02,  -9.95573770e+02],
       [ -1.00000000e+03,  -9.96721311e+02,   1.14754098e+00,
         -9.95573770e+02,  -9.91147541e+02]])

Secondly, the matrix in your numpy case is practically singular, its determinant is 1e-21. This means that its inverse doesn't exist, and is horribly ill-defined at best (explaining the matrix elements of magnitude 1e16) in the result.

Note that the difference between the matrices cannot be handwaved away saying that "internally the two matrices are the same, only matlab didn't write as many decimal places", since a matrix element of numpy 0.13535533905932737 is printed as 0.1340 in matlab, which is clearly wrong, irrespective of printing precision.

Also note that the only reason your matlab matrix is not that singular is that most of its elements are 0.135, but its first element is 0.134. If you make this element become equal to the others (just like in the numpy case!), you get

>> b=a
>> b(1,1)=0.135
>> det(b)
ans =
     1.481453848484194e-21

There you have it. If your actual matrix is the one in numpy, then the solution is easy: don't invert it, since it doesn't have an inverse.

---

To convince you, consider a random matrix with the same pattern of elements as your newK matrix:

k1=np.random.rand()*2-1 #element [0,0]
k2=np.random.rand()*2-1 #element [0,1]
k3=k1+k2 #element [0,-1]

Krand=np.array([[k1,k2,0,0,k3],[k2,k1,0,k3,0],[0,0,k1,-k2,k2],[0,k3,-k2,k1,k2],[k3,0,k2,k2,k1]])

then you get

In [322]: np.linalg.det(Krand)
Out[322]: 3.3121334687895703e-17

meaning that any matrix with this structure will be singular.

Final proof using sympy:

import sympy as sym
k1,k2=sym.symbols('k1,k2')
Ksym=sym.Matrix([[k1,k2,0,0,k1+k2],[k2,k1,0,k1+k2,0],[0,0,k1,-k2,k2],[0,k1+k2,-k2,k1,k2],[k1+k2,0,k2,k2,k1]])

Then the symbolic determinant is

In [347]: sym.det(Ksym)
Out[347]: 0