# Chapter 7 - Direct solutions of linear equations¶

## Example No. 7_01 Pg No. 211¶

In [1]:
from numpy import mat
#Elimination Process

A = mat([[3, 2, 1, 10],[2, 3, 2, 14],[1, 2, 3, 14]])
A[1,:] = A[1,:] - A[0,:]*A[1,0]/A[0,0]
A[2,:] = A[2,:] - A[0,:]*A[2,0]/A[0,0]
print 'A =\n',A

print 'After transformation'
A[2,:] = A[2,:] - A[1,:]*A[2,1]/A[1,1]
print 'A=\n',A

z = A[2,3]/A[2,2]
y = (A[1,3] - A[1,2]*z)/A[1,1]
x = (A[0,3] - A[0,1]*y - A[0,2]*z)/A[0,0]
print 'x =',x
print 'y =',y
print 'z =',z

A =
[[ 3  2  1 10]
[ 0  2  2  8]
[ 0  2  3 11]]
After transformation
A=
[[ 3  2  1 10]
[ 0  2  2  8]
[ 0  0  1  3]]
x = 1
y = 1
z = 3


## Example No. 7_02 Pg No. 214¶

In [2]:
from numpy import mat,shape
#Basic Gauss Elimination

A = mat([[  3,  6,  1],[ 2,  4 , 3],[ 1,  3,  2 ]])
B = [16, 13, 9]
S=shape(A)
ar1,ac1 = S[0],S[1]
Aug = mat([[3, 6 ,1 ,16],[2, 4, 3, 13],[1, 3, 2, 9]])
for i in range(1,ar1):
Aug[i,:] = Aug[i,:] - (Aug[i,0]/Aug[0,0])*Aug[0,:]

print Aug
print 'since Aug(2,2) = 0 elimination is not possible,so reordering the matrix'
temp=A[2,:]
A[2,:]=A[1,:]
A[1,:]=temp
print Aug
print 'Elimination is complete and by back substitution the solution is\n'
print 'x3 = 1, x2 = 2 , x1 = 1 '

[[ 3  6  1 16]
[ 2  4  3 13]
[ 1  3  2  9]]
since Aug(2,2) = 0 elimination is not possible,so reordering the matrix
[[ 3  6  1 16]
[ 2  4  3 13]
[ 1  3  2  9]]
Elimination is complete and by back substitution the solution is

x3 = 1, x2 = 2 , x1 = 1


## Example No. 7_03 Pg No. 220¶

In [1]:
from numpy import array , zeros, dot, diag

A = array([[ 2.,  2.,  1.],[4.,  2.,  3.],[ 1.,  -1.,  1.]])
B = array([[6.],[4.],[0.]])

## Gauss Elimination using partial pivoting

def GEPP(A, b):
'''
Gaussian elimination with partial pivoting.
% input: A is an n x n nonsingular matrix
%        b is an n x 1 vector
% output: x is the solution of Ax=b.
% post-condition: A and b have been modified.
'''
n =  len(A)
if b.size != n:
raise ValueError("Invalid argument: incompatible sizes between A & b.", b.size, n)
# k represents the current pivot row. Since GE traverses the matrix in the upper
# right triangle, we also use k for indicating the k-th diagonal column index.
for k in xrange(n-1):
#Choose largest pivot element below (and including) k
maxindex = abs(A[k:,k]).argmax() + k
if A[maxindex, k] == 0:
raise ValueError("Matrix is singular.")
#Swap rows
if maxindex != k:
A[[k,maxindex]] = A[[maxindex, k]]
b[[k,maxindex]] = b[[maxindex, k]]
for row in xrange(k+1, n):
multiplier = A[row][k]/A[k][k]
#the only one in this column since the rest are zero
A[row][k] = multiplier
for col in xrange(k + 1, n):
A[row][col] = A[row][col] - multiplier*A[k][col]
#Equation solution column
b[row] = b[row] - multiplier*b[k]
#print A ;print b
x = zeros(n)
k = n-1
x[k] = b[k]/A[k,k]
while k >= 0:
x[k] = (b[k] - dot(A[k,k+1:],x[k+1:]))/A[k,k]
k = k-1
return x
Aug=GEPP(A,B)
print 'Gauss Elimination using partial pivoting : ans is\n',Aug

#Back Substitution
def forward_elimination(A, b, n):
"""
Calculates the forward part of Gaussian elimination.
"""
for row in range(0, n-1):
for i in range(row+1, n):
factor = A[i,row] / A[row,row]
for j in range(row, n):
A[i,j] = A[i,j] - factor * A[row,j]

b[i] = b[i] - factor * b[row]

return A, b

def back_substitution(a, b, n):
""""
Does back substitution, returns the Gauss result.
"""
x = zeros((n,1))
x[n-1] = b[n-1] / a[n-1, n-1]
for row in range(n-2, -1, -1):
sums = b[row]
for j in range(row+1, n):
sums = sums - a[row,j] * x[j]
x[row] = sums / a[row,row]
return x

def gauss(A, b):
"""
This function performs Gauss elimination without pivoting.
"""
n = A.shape[0]

# Check for zero diagonal elements
if any(diag(A)==0):
raise ZeroDivisionError(('Division by zero will occur; '
'pivoting currently not supported'))

A, b = forward_elimination(A, b, n)
return back_substitution(A, b, n)

x = gauss(A, B)
print ('\n The solution can be obtained by back substitution \n x1 = %i \n x2 = %i \n x3 = %i \n'%(x[0], x[1], x[2]))

Gauss Elimination using partial pivoting : ans is
[  9.  -1. -10.]

The solution can be obtained by back substitution
x1 = 4
x2 = 0
x3 = -4



## Example No. 7_04 Pg No. 228¶

In [1]:
from numpy import mat, shape
#Gauss Jordan Elimination

A = mat([[  2,  4,  -6],[1,  3,  1],[2,  -4,  -2  ]])
B = mat([[  -8],[ 10],[ -12  ]])
S=shape(A)
ar, ac =S[0],S[1]
Aug = mat([[  2,  4,  -6,  -8],[ 1,  3,  1,  10],[ 2,  -4,  -2,  -12  ]])
print 'Aug = \n',Aug

for i in range(0,ar):
Aug[i,i:ar+1] = Aug[i,i:ar+1]/Aug[i,i]
print 'Aug = \n',Aug
for k in range(0,ar):
if k != i:
Aug[k,i:ar] = Aug[k,i:ar] - Aug[k,i]*Aug[i,i:ar]

print 'Aug = \n',Aug

Aug =
[[  2   4  -6  -8]
[  1   3   1  10]
[  2  -4  -2 -12]]
Aug =
[[  1   2  -3  -4]
[  1   3   1  10]
[  2  -4  -2 -12]]
Aug =
[[  1   2  -3  -4]
[  0   1   4  10]
[  0  -8   4 -12]]
Aug =
[[  1   2  -3  -4]
[  0   1   4  10]
[  0  -8   4 -12]]
Aug =
[[  1   0 -11  -4]
[  0   1   4  10]
[  0   0  36 -12]]
Aug =
[[  1   0 -11  -4]
[  0   1   4  10]
[  0   0   1  -1]]
Aug =
[[ 1  0  0 -4]
[ 0  1  0 10]
[ 0  0  1 -1]]


## Example No. 7_05 Pg No. 234¶

In [2]:
from numpy import array, arange, zeros
from scipy.linalg import lu
from __future__ import division

#DoLittle LU Decomposition

A = array([[ 3,  2,  1],[2 , 3 , 2],[1,  2,  3  ]])
B = array([[ 10],[14],[14  ]])
n,n = A.shape
L=lu(A)[1]
U=lu(A)[2]
print 'U = \n',U
print 'L = \n',L

z=zeros([3,1])
#Forward Substitution
for i in range(0,n):
z[i,0] = B[i,0]
for j in range(0,i-1):
z[i,0] = z[i,0] - L[i,0]*z[j,0];

print '\n z(%i) = %f \n'%(i+1,z[i,0])

x=zeros([3,1])
#Back Substitution
for i in arange(n-1,-1,-1):
x[i,0] = z[i,0]
for j in arange(n-1,i+1,-1):
x[i,0] = x[i,0] - U[i,j]*x[j,0]

x[i,0] = x[i,0]/U[i,i]
print '\n x(%i) = %f \n'%(i+1,x[i,0])

U =
[[ 3.          2.          1.        ]
[ 0.          1.66666667  1.33333333]
[ 0.          0.          1.6       ]]
L =
[[ 1.          0.          0.        ]
[ 0.66666667  1.          0.        ]
[ 0.33333333  0.8         1.        ]]

z(1) = 10.000000

z(2) = 14.000000

z(3) = 10.666667

x(3) = 6.666667

x(2) = 8.400000

x(1) = 1.111111



## Example No. 7_06 Pg No. 242¶

In [1]:
from numpy import sqrt, mat,shape, zeros
#Cholesky's Factorisation

A = mat([[ 1, 2, 3],[2, 8 ,22],[3, 22, 82 ]])
n= shape(A)[0]
U=zeros([n,n])

for i in range(0,n):
for j in range(0,n):
if i == j:
U[i,i] = A[i,i]
for k in range(0,i-1):
U[i,i] = U[i,i]-U[k,i]**2

U[i,i] = sqrt(U[i,i])
elif i < j:
U[i,j] = A[i,j]
for k in range(0,i-1):
U[i,j] = U[i,j] - U[k,i]*U[k,j]

U[i,j] = U[i,j]/U[i,i]

print 'U =\n',U

U =
[[ 1.          2.          3.        ]
[ 0.          2.82842712  7.77817459]
[ 0.          0.          8.54400375]]


## Example No. 7_07 Pg No. 245¶

In [2]:
from numpy import mat
#Ill-Conditioned Systems

A = mat([[ 2,  1],[2.001,  1]])
B = mat([[ 25],[25.01 ]])
x=[0,0]
x[0] = (25 - 25.01)/(2 - 2.001)
x[1] = (25.01*2 - 25*2.001)/(2*1 - 2.001*1)
print '\n x(1) = %f \n x(2) = %f \n'%(x[1],x[1])
x[0] = (25 - 25.01)/(2 - 2.0005)#
x[1] = (25.01*2 - 25*2.0005)/(2*1 - 2.0005*1)#
print '\n x(1) = %f \n x(2) = %f \n'%(x[1],x[1])
x=mat([[x[0]],[x[1]]])
r = A*(x)-B
print 'x=\n',x
print 'r=\n',r

 x(1) = 5.000000
x(2) = 5.000000

x(1) = -15.000000
x(2) = -15.000000

x=
[[ 20.]
[-15.]]
r=
[[ -2.66453526e-12]
[  1.00000000e-02]]