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
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 '
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]))
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
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])
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
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