Chapter 10 : LU Decomposition and matrix inverse

Ex10.1 Page 277

In [9]:
from numpy import mat
from numpy.linalg import det
A = mat([[3,-0.1,-0.2],[0.1,7,-0.3],[0.3,-0.2,10]])
U = A#
print "A =",A
m = U[0,0]
n = U[1,0]
a = n/m#
U[1:2] = U[1:2] - U[0:1] / (m/n)#
n = U[2,0]
b = n/m

U[2:3] = U[2:3] - U[0:1] / (m/n)#
m = U[1,1]
n = U[2,1]
c = n/m#
U[2:3] = U[2:3] - U[1:2] / (m/n)#
print "\nU = \n",U
L = mat([[1,0,0],[a,1,0],[b,c,1]])
print "\nL calculated based on gauss elimination method = \n",L

A = [[  3.   -0.1  -0.2]
 [  0.1   7.   -0.3]
 [  0.3  -0.2  10. ]]

U = 
[[  3.          -0.1         -0.2       ]
 [  0.           7.00333333  -0.29333333]
 [  0.           0.          10.01204188]]

L calculated based on gauss elimination method = 
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]

Ex10.2 Page 279

In [11]:
from numpy import mat
A = mat([[3,-0.1,-0.2],[0.1,7,-0.3],[0.3,-0.2,10]])
B = mat([[7.85],[-19.3],[71.4]])
X = (A**-1) * B
print "X = \n",X

X = 
[[ 3. ]
 [-2.5]
 [ 7. ]]

Ex10.3 Page 284

In [24]:
from numpy import mat,array
A = mat([[3,-0.1,-0.2],[0.1,7,-0.3],[0.3,-0.2,10]])
B = (A**-1)
L = mat([[1,0,0],[0.033333,1,0],[0.1,-0.02713,1]])
U = mat([[3,-0.1,-0.2],[0,7.0033,-0.293333],[0,0,10.012]])
for i  in range(1,4):
    if i==1:
        m = mat([[1],[0],[0]])
    else:
        if i==2:
            m = mat([[0],[1],[0]])
        else:
            m = mat([[0],[0],[1]])
        
    
    d = (L**-1) * m#
    x = (U**-1) * d#
    B[:,i-1] = x

print "B=\n",B

B=
[[ 0.33248872  0.00494409  0.00679813]
 [-0.00518174  0.14290333  0.00418348]
 [-0.01007834  0.00270975  0.09988014]]

Ex10.4 Page 291

In [32]:
from numpy import mat
from __future__ import division
A = mat([[1,1/2,1/3],[1/2,1/3,1/4],[1/3,1/4,1/5]])
n = A[1,0]
A[1,:] = A[1,:]/n
n = A[2,0]
A[2,:] = A[2,:]/n
B = (A**-1)#
print "A = \n",A

m=range(1,4)
su=range(1,4)
for j in range(1,4):
    a = 0#
    for i in range(1,4):
        m[i-1]= A[j-1,i-1]
        su[j-1] = a + m[i-1]#
        a = su[j-1]#


if su[0]< su[1]:
    if su[1]< su[2]:
        z = su[2]
    else:
        z = su[1]
    
else:
     if su[0] < su[2]:
        z = su[2]#
     else:
        z = su[0]#
    
m=range(1,4)
sm=range(1,4)
for j in  range(1,4):
    a = 0#
    for i in range(1,4):
        m[i-1]= B[j-1,i-1]
        sm[j-1]= a + abs(m[i-1])
        a = sm[j-1]#


if sm[0]< sm[1]:
    if sm[1]< sm[2] :
        y = sm[2]
    else:
        y = sm[1]
    
else:
     if sm[0]< sm[2]:
        y = sm[2]
     else:
        y = sm[0]#
    

C = z*y#
print "\nCondition number for the matrix =\n",C

A = 
[[ 1.          0.5         0.33333333]
 [ 1.          0.66666667  0.5       ]
 [ 1.          0.75        0.6       ]]

Condition number for the matrix =
451.2