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 [ ]: