Numeric Analysis 4 - Numeric Linear Algebra

Numeric Linear Algebra¶

---

Gauss Elimination¶

[Example] Gauss Elimination.

In [1]:

def Gauss( A ):
    n = len(A)     # number of rows
    m = len(A[0])  # number of columns
     
    for i in range(n):
        
        # Pivoting
        
        # find pivot row   pivot =the largest absolute value among A[k][i]
        Amax = 0.0
        kmax = i
        for k in range(i+1,n):
            if abs(A[k][i]) > Amax :
                Amax = abs(A[k][i])
                kmax = k
        
        # swap row i and row kmax
        if kmax != i : A[i], A[kmax] = A[kmax], A[i]
        
        
        # Row operations
        
        aa = A[i][i]
        # make diagonal element unity dividing the row by the diagonal element
        for j in range(i,m): A[i][j] /= aa
        
        # make off-diagonal element zero below row i
        for k in range(i+1,n):
            aa = A[k][i]
            for j in range(i,m):
                A[k][j] += A[i][j]*(-aa)
        
    
    # back substitution
    x = [0]*n
    
    for i in range(n-1,-1,-1):   # i = n-1, ... , 1, 0
        x[i] = A[i][-1]          # last column
        for j in range(i+1,n):
            x[i] -= A[i][j] * x[j]

    return x

In [2]:

A = [ [0,8,2,-7], 
      [3,5,2,8], 
      [6,2,8,26] ]

sol = Gauss( A )

n = len(A)      
m = len(A[0]) 
for i in range(n):
    for j in range(m):
        print('%8.4f ' % A[i][j], end='' )
    print()

print('\n solution = ', sol )

  1.0000   0.3333   1.3333   4.3333 
  0.0000   1.0000   0.2500  -0.8750 
  0.0000   0.0000   1.0000   0.5000 

 solution =  [3.9999999999999996, -1.0, 0.5]

Gauss–Jordan Elimination¶

The Inverse of a Matrix¶

[Example] from sec. 7.8

In [3]:

def Gauss_Jordan( A ):
    n = len(A)     # number of rows
    m = len(A[0])  # number of columns
     
    for i in range(n):
        
        # Pivoting
        
        # find pivot row   pivot =the largest absolute value among A[k][i]
        Amax = 0.0
        kmax = i
        for k in range(i+1,n):
            if abs(A[k][i]) > Amax :
                Amax = abs(A[k][i])
                kmax = k
        
        # swap row i and row kmax
        if kmax != i : A[i], A[kmax] = A[kmax], A[i]
        
        
        # Row operations
        
        aa = A[i][i]
        # make diagonal element unity dividing the row by the diagonal element
        for j in range(i,m): A[i][j] /= aa
        
        # make off-diagonal element zero above and below row i
        for k in range(n):
            if k == i : continue
            aa = A[k][i]
            for j in range(i,m):
                A[k][j] += A[i][j]*(-aa)
            
    return 

In [4]:

A = [ [-1,1,2], 
      [3,-1,1], 
      [-1,3,4] ]

n = len(A) 

# pad unit matrix on the right
for i in range(n):
    for j in range(n):
        if i==j : 
            A[i].append(1)
        else :
            A[i].append(0)

m = len(A[0])
for i in range(n):
    for j in range(m):
        print('%7.3f ' % A[i][j], end='' )
    print()

    
Gauss_Jordan( A )


print('\n inverse matrix is \n')
for i in range(n):
    for j in range(m):
        print('%7.3f ' % A[i][j], end='' )
    print()     

 -1.000   1.000   2.000   1.000   0.000   0.000 
  3.000  -1.000   1.000   0.000   1.000   0.000 
 -1.000   3.000   4.000   0.000   0.000   1.000 

 inverse matrix is 

  1.000   0.000   0.000  -0.700   0.200   0.300 
  0.000   1.000   0.000  -1.300  -0.200   0.700 
  0.000   0.000   1.000   0.800   0.200  -0.200