Solutions for Lab Class on Linear Systems and Conditioning

In [6]:
# 1(a)
# Code for testing the steps of the Gauss-Jordan method
# for finding the solution 'x' to a linear system A x = b,
# along with the inverse matrix A^{-1}.

import numpy as np
np.set_printoptions(threshold=np.inf) # print the whole augmented array

def swaprows(M,i,j):
    M[i-1,:], M[j-1,:] = M[j-1,:].copy(), M[i-1,:].copy()

def augmented(A,b):
    n = A.shape[0]
    aug = np.zeros((n,2*n+1))
    aug[:,:n] = A
    aug[:,n:2*n] = np.identity(n)
    aug[:,2*n] = b.transpose()
    return aug

A = np.matrix("0,2,1;2,-1,1;1,3,2", dtype=np.float64)
b = np.matrix("7;3;13", dtype=np.float64)
M = augmented(A, b)
print(M,"\n")  # Show the augmented matrix

swaprows(M,1,2) # Swap the first and second rows
print(M,"\n")

M[0,:] = M[0,:] / 2  # Divide the first row by 2
print(M,"\n")

M[2,:] = M[2,:] - M[0,:]  # Subtract the 1st row from the 3rd row
print(M,"\n")

M[2,:] = M[2,:]*4
M[2,:] = M[2,:] - 7*M[1,:]
M[1,:] = M[1,:]/2
print(M,"\n")

M[1,:] = M[1,:] + M[2,:]/2
M[0,:] = M[0,:] + M[2,:]/2
M[2,:] = -M[2,:]
print(M,"\n")

M[0,:] = M[0,:] + M[1,:]/2
print(M,"\n")

print("The matrix inverse is\n", M[:,3:6], "\nand the solution is\n x = ", M[:,-1])
[[  0.   2.   1.   1.   0.   0.   7.]
 [  2.  -1.   1.   0.   1.   0.   3.]
 [  1.   3.   2.   0.   0.   1.  13.]] 

[[  2.  -1.   1.   0.   1.   0.   3.]
 [  0.   2.   1.   1.   0.   0.   7.]
 [  1.   3.   2.   0.   0.   1.  13.]] 

[[  1.   -0.5   0.5   0.    0.5   0.    1.5]
 [  0.    2.    1.    1.    0.    0.    7. ]
 [  1.    3.    2.    0.    0.    1.   13. ]] 

[[  1.   -0.5   0.5   0.    0.5   0.    1.5]
 [  0.    2.    1.    1.    0.    0.    7. ]
 [  0.    3.5   1.5   0.   -0.5   1.   11.5]] 

[[ 1.  -0.5  0.5  0.   0.5  0.   1.5]
 [ 0.   1.   0.5  0.5  0.   0.   3.5]
 [ 0.   0.  -1.  -7.  -2.   4.  -3. ]] 

[[ 1.  -0.5  0.  -3.5 -0.5  2.   0. ]
 [ 0.   1.   0.  -3.  -1.   2.   2. ]
 [-0.  -0.   1.   7.   2.  -4.   3. ]] 

[[ 1.  0.  0. -5. -1.  3.  1.]
 [ 0.  1.  0. -3. -1.  2.  2.]
 [-0. -0.  1.  7.  2. -4.  3.]] 

The matrix inverse is
 [[-5. -1.  3.]
 [-3. -1.  2.]
 [ 7.  2. -4.]] 
and the solution is
 x =  [ 1.  2.  3.]
In [13]:
# 1(b)
A = np.matrix("0,2,3;0,1,2;1,0,5", dtype=np.float64)
b = np.matrix("1;1;7", dtype=np.float64)
M = augmented(A, b)

# Hint: start by swapping rows 1 and 3. 
swaprows(M,1,3)
print(M,"\n")

M[2,:] =  M[2,:] - 2*M[1,:]
M[2,:] = -M[2,:]
print(M,"\n")

M[1,:] = M[1,:] - 2*M[2,:]
M[0,:] = M[0,:] - 5*M[2,:]
print(M,"\n")

print("The matrix inverse is\n", M[:,3:6], "\nand the solution is\n x = ", M[:,-1])
[[ 1.  0.  5.  0.  0.  1.  7.]
 [ 0.  1.  2.  0.  1.  0.  1.]
 [ 0.  2.  3.  1.  0.  0.  1.]] 

[[ 1.  0.  5.  0.  0.  1.  7.]
 [ 0.  1.  2.  0.  1.  0.  1.]
 [-0. -0.  1. -1.  2. -0.  1.]] 

[[  1.   0.   0.   5. -10.   1.   2.]
 [  0.   1.   0.   2.  -3.   0.  -1.]
 [ -0.  -0.   1.  -1.   2.  -0.   1.]] 

The matrix inverse is
 [[  5. -10.   1.]
 [  2.  -3.   0.]
 [ -1.   2.  -0.]] 
and the solution is
 x =  [ 2. -1.  1.]

1(c). See gauss-jordan3.py.

Question 3

In [21]:
# 3(a)
def row_sum_norm(A):
    rows,cols = A.shape
    B = np.abs(A)
    rowsums = [B[i,:].sum() for i in range(rows)]
    return np.max(rowsums)

def condition_number(A):
    t0 = row_sum_norm(A)
    t1 = row_sum_norm(np.linalg.inv(A))
    return t0*t1
In [26]:
# 3(b)
A1 = np.array([[1,2],[3,4]])
A2 = np.array([[1,1],[1,1.001]])
A3 = np.array([[1,1/2,1/3],[1/2,1/3,1/4],[1/3,1/4,1/5]])

print("Condition numbers:")
print(condition_number(A1))
print(condition_number(A2))
print(condition_number(A3))

# Alternatively, these can be calculated using numpy. Note the optional parameter p = 1 here.
print(np.linalg.cond(A1, p=1))
print(np.linalg.cond(A2, p=1))
print(np.linalg.cond(A3, p=1))
Condition numbers:
21.0
4004.001
748.0
21.0
4004.001
748.0
In [ ]: