Gaussian Elimination to solve linear and non-linear system of equations.

gausse.py

"""
Gaussian Elimination with Partial Pivoting.

This module contains 5 functions which procedurally solve
the linear system Ax = b, A is an nxn matrix and b is a
column vector with n rows.
"""
import numpy
import math


 
def pivot(A):
    """
    Pivots matrix A - finds row with maximum first entry
    and if nessecary swaps it with the first row.
    Input Arguments
    ---------------
    Augmented Matrix A

    Returns
    -------
    Pivoted Augmented Matrix A
    """
    
    B = numpy.zeros((1,2))
    B[0,0]=A.shape[0]       
    B[0,1]=A.shape[1]
    nrows =B[0,0] #This stores dimensions of the 
    ncols =B[0,1] #matrix in an array
    
    pivot_size = numpy.abs(A[0,0])        
    #checks for largest first entry and swaps
    pivot_row = 0; 
    for i in range(int(0),int(nrows)):
        if numpy.abs(A[i,0])>pivot_size:
            pivot_size=numpy.abs(A[i,0])
            pivot_row = i
        
    if pivot_row>0:
        tmp = numpy.empty(ncols)
        tmp[:] = A[0,:]
        A[0,:] = A[pivot_row,:]
        A[pivot_row,:] = tmp[:]
        
    return A
# Testing the function
if __name__ == "__main__": 

    A = numpy.array([[-1.,1.,2.,1.],[3.,-1.,1.,1.],[-1.,3.,4.,1.]])
    print 'Testing pivot'
    print 'A before pivot\n',A
    print 'A after pivot\n',pivot(A),'\n'


def elim(A):
    """
    elim(A)uses row operations to introduce
    zeros below the diagonal in first column
    of matrix A
    Input Argument
    ---------------
    Augmented Matrix A

    Returns
    -------
    A with eliminated first column
    """
    A = pivot(A)
    B = numpy.zeros((1,2))
    B[0,0]=A.shape[0]       
    B[0,1]=A.shape[1]
    nrows =B[0,0]
    ncols =B[0,1]
    
    #row operations
    rpiv = 1./float(A[0][0])
    
    for irow in range(int(1),int(nrows)):
        s=A[irow,0]*rpiv
        for jcol in range(int(0),int(ncols)):
            A[irow,jcol] = A[irow,jcol] - s*A[0,jcol]
        
    return A
# Testing the function
if __name__ == "__main__": 

    A = numpy.array([[-1.,1.,2.,1.],[3.,-1.,1.,1.],[-1.,3.,4.,1.]])
    print 'Testing elimination'
    print 'A before elimination\n',pivot(A)
    print 'A after elimination\n',elim(pivot(A)),'\n'



def backsub(A):
    """
    backsub(A) solves the upper triangular system

    Input Argument
    ---------------
    Augmented Matrix A

    Returns
    -------
    vector b, solution to the linear system
    """
    B = numpy.zeros((1,2))
    B[0,0]=A.shape[0]       
    B[0,1]=A.shape[1]
    n =B[0,0]
    ncols =B[0,1]
        
    x=numpy.zeros((n,1))

    x[n-1]=A[n-1,n]/A[n-1,n-1]
    for i in range(int(n-1),int(-1),int(-1)):
        for j in range(int(i+1),int(n)):
            A[i,n] = A[i,n]-A[i,j]*x[j]

        x[i] = A[i,n]/A[i,i]
    return x

# Testing the function
if __name__ == "__main__": 

    A = numpy.array([[-1.,1.,2.,1.],[3.,-1.,1.,1.],[-1.,3.,4.,1.]])
    print 'Testing backsub'
    print 'A before backsub\n',A,'\n'
    print 'A after backsub\n',backsub(A),'\n'




def gaussfe(A):
    """
    gaussfe(A)uses row operations to reduce A
    to upper triangular form by calling elim and
    pivot
    Input Argument
    ---------------
    Augmented Matrix A

    Returns
    -------
    A in upper triangular form
    """
    B = numpy.zeros((1,2))
    B[0,0]=A.shape[0]       
    B[0,1]=A.shape[1]
    nrows =B[0,0]
    ncols =B[0,1]

    for i in range(int(0),int(nrows-1)):
        A[i:nrows,i:ncols]=pivot(A[i:nrows,i:ncols])
        A[i:nrows,i:ncols]=elim(A[i:nrows,i:ncols])

    return A
# Testing the function
if __name__ == "__main__": 

    A = numpy.array([[-1.,1.,2.,1.],[3.,-1.,1.,1.],[-1.,3.,4.,1.]])
    print 'Testing gaussfe'
    print 'A before gaussian elimination\n',gaussfe(pivot(A))
    print 'A after g.elimination\n',gaussfe(elim(pivot(A))),'\n'
    

def solve(A,b):
    """
    Solve augments the nxn matrix A and column vector b
    then calls upon the functions in this module to solve
    the linear system

    Input Argument
    ---------------
    nxn matrix A and column vector b

    Returns
    -------
    x, the solution to the linear system
    """
    B = numpy.zeros((1,2))
    B[0,0]=A.shape[0]
    B[0,1]=A.shape[1]
    nrows =B[0,0]
    ncols =B[0,1]
    Aug= numpy.zeros((nrows,ncols+1))
    #these 2 loops augment the matrix with column vector
    for i in range(int(0),int(nrows)):
        for j in range(int(0),int(ncols)):
            Aug[i,j] = A[i,j]

    for k in range(int(0),int(nrows)):
        
        Aug[k,ncols]= b[k]

    A = Aug

    A=gaussfe(A)
    x=backsub(A)
    x=x.T
    return x

# Specific test
if __name__ == "__main__": 
    A = numpy.array([[1.,1.,1.],[1.,-2.,2.],[1.,2.,-1.]])
    b = numpy.array([[0.],[4.],[2.]])
    print 'Testing solve'
    print 'A is \n',A
    print 'B is \n',b
    print 'The solution is \n', solve(A,b),'\n\n\n'


    

# Testing the function generally
if __name__ == "__main__": 
    import numpy
    import time
    print "Tests for solve\n\n"
    
    n = 100
    print "Matrix size",n
    A = numpy.random.random((n,n))
    xe = numpy.random.random(n)
    b = numpy.dot(A,xe)
    solve_start = time.clock()
    x=solve(A,b)
    print "solve time ",str(time.clock()-solve_start)
    err = numpy.max(numpy.abs(xe-x))
    print "Error is ",err,"\n\n"
    assert err<1e-10    
    print "solve test passed."

newton.py

"""
newton computes the jacobian matrix of a non linear system
and uses newton method to solve the non linear system


"""
import numpy
import math
import gausse

def jacobian(f,x):
    """
    jacobian works out jacobian matrix of non linear system
    
    Input Arguments
    ---------------
    Function containing non linear system f, initial guess x
    
    Returns
    ------
    Jacobian values
    """
    h=1.0e-4
    n=len(x)
    jac=numpy.zeros((n,n))
    temp = numpy.zeros((n,1))
    f0=f(x) # using finite difference to calculate jacobian
    for i in range(n):
        temp[:]=x[:]
        temp[i] = temp[i]+h
        f1=f(temp)
        jac[:,i] = (f1-f0)/h
    return jac, f0



def newtonraphson(f,x,tol=1.0e-9):
    """
    newtonraphson uses newtons method to solve non linear systems

    Input Arguments
    ---------------
    Vector of functions f, initial guess values x, error tollerance

    Returns
    -------
    Solutions to non linear system
    """
    
    for i in range(50):
               
        jac,f0=jacobian(f,x)
        print "jac and f0\n\n",jac,"\n\n", f0,"\n\n"
        dx = gausse.solve(jac,f0)
        print 'dx is',(dx)
        dx=dx.T
        x=x-dx
        print 'x is', x
        print "max is  ",max(dx)
        if max(abs(dx)) < tol: return x
    print ' too many interations'





#Tessting the function#

if __name__ == "__main__":
    def f(x):
        f = numpy.zeros(len(x))
        f[0] = x[0]**2 - 2
        f[1] = x[1]**2 - 2
   
        return f

    x= numpy.zeros((2,1))
    x[0]=1
    x[1]=1
    print 'with the initial value as 1,1 we get \n',newtonraphson(f,x,tol=1.0e-9),'\n'
    NLsol=newtonraphson(f,x,tol=1.0e-9)
    for k in range(0,2):
        if 1.4<numpy.abs(NLsol[k])<1.42: #the approximate bound for the solution
            print 'test passed for root',k,'\n'
            
    x= numpy.zeros((2,1))
    x[0]=-1
    x[1]=-1
    print 'with the initial value as -1,-1 we get \n',newtonraphson(f,x,tol=1.0e-9),'\n'
    NLsol=newtonraphson(f,x,tol=1.0e-9)
    for k in range(0,2):
        if 1.4<numpy.abs(NLsol[k])<1.42: #the approximate bound for the solution
            print 'test passed for root',k,'\n'
    
    x= numpy.zeros((2,1))
    x[0]=1
    x[1]=-1
    print 'with the initial value as 1,-1 we get \n',newtonraphson(f,x,tol=1.0e-9),'\n'
    NLsol=newtonraphson(f,x,tol=1.0e-9)
    for k in range(0,2):
        if 1.4<numpy.abs(NLsol[k])<1.42: #the approximate bound for the solution
            print 'test passed for root',k,'\n'

    x= numpy.zeros((2,1))
    x[0]=-1
    x[1]=1
    print 'with the initial value as -1,1 we get \n',newtonraphson(f,x,tol=1.0e-9),'\n'
    NLsol=newtonraphson(f,x,tol=1.0e-9)
    for k in range(0,2):
        if 1.4<numpy.abs(NLsol[k])<1.42: #the approximate bound for the solution
            print 'test passed for root',k,'\n'

q3system.py

import numpy
import math
import newton

#The first system of non-linear equations
def f(x):
    f = numpy.zeros(len(x))
    f[0] = x[0]**2 + x[1]**2 - 1
    f[1] = 5.0*x[0]**2 + 21*x[1]**2 - 9
   
    return f
x= numpy.zeros((2,1))
x[0]=0.5
x[1]=0.5
print "Solutions are \n",newton.newtonraphson(f,x,tol=1.0e-9),"\n"

def f(x):
    f = numpy.zeros(len(x))
    f[0] = x[0]**2 + x[1]**2 - 1
    f[1] = 5.0*x[0]**2 + 21*x[1]**2 - 9
   
    return f
x= numpy.zeros((2,1))
x[0]=-0.5
x[1]=-0.5
print newton.newtonraphson(f,x,tol=1.0e-9),"\n"

def f(x):
    f = numpy.zeros(len(x))
    f[0] = x[0]**2 + x[1]**2 - 1
    f[1] = 5.0*x[0]**2 + 21*x[1]**2 - 9
   
    return f
x= numpy.zeros((2,1))
x[0]=-0.5
x[1]=0.5
print newton.newtonraphson(f,x,tol=1.0e-9),"\n"

def f(x):
    f = numpy.zeros(len(x))
    f[0] = x[0]**2 + x[1]**2 - 1
    f[1] = 5.0*x[0]**2 + 21*x[1]**2 - 9
   
    return f
x= numpy.zeros((2,1))
x[0]=0.5
x[1]=-0.5
print newton.newtonraphson(f,x,tol=1.0e-9),"\n"




#The second system of non-linear equations
def f(x):
    
    f = numpy.zeros(len(x))
    f[0] = x[0]**2 + x[1]**2 + x[2]**2 - 1
    f[1] = 2*x[0]**2 + x[1]**2 - 4*x[2]
    f[2] = 3*x[0]**2 - 4*x[1] + x[2]**2
   
    return f
x= numpy.zeros((3,1))
x[0]=1
x[1]=1
x[2]=1
print "Solution is \n",newton.newtonraphson(f,x,tol=1.0e-9)