Solving a 9x9 matric with gausian emlinination with pivoting in python

Question

Im needing to solve a whole range of 8x8 and 9x9 matrices so thought I could build a python program to make the whole thing easier.

So far I have managed to create:

from __future__ import division
import numpy as np

def solveEqns(A,v):
    def lu( A ):
        #Factor A into LU by Gaussian elimination with scaled partial pivoting
        n, m = np.shape( A )
        if n != m:
            print "Error: input matrix is not square"
            return None
            # Generate initial index vector
        p = range( n )
            # Determine the largest (in magnitude) element in each row.  These
            # factors are used to scale the pivot elements for comparison purposes
            # when deciding which row to use as a pivot row.
        s = [0] * n
        for i in xrange( n ):
            smax = 0.0
            for j in xrange( n ):
                smax = max( smax, abs( A[i][j] ) )
            s[i] = smax
            # Begin Gaussian elimination.
        for k in xrange( n - 1 ):
        # Find the remaining row with the largest scaled pivot.
            rmax = 0.0
            for i in xrange( k, n ):
                r = abs( A[p[i][k]] / s[p[i]] )
                if r > rmax:
                    rmax = r
                    j = i
        # Row j has the largest scaled pivot, so "swap" that row with the
        # current row (row k).  The swap is not actually done by copying rows,
        # but by swaping two entries in an index vector.
            p[j], p[k] = ( p[k], p[j] )
        # Now carry out the next elimination step as usual, except for the
        # added complication of the index vector.
            for i in xrange( k + 1, n ):
                xmult = A[p[i],k] / A[p[k],k]
                A[p[i],k] = xmult
                for j in xrange( k + 1, n ):
                    A[p[i],j] = A[p[i],j] - xmult * A[p[k],j]
    # All done, return factored matrix A and permutation vector p
        return ( A, p )

def solve( A, p, b ):
#Solves Ax = b given an LU factored matrix A and permuation vector p
    n, m = np.shape( A )
    if n != m:
        print "Error: input matrix is not square"
        return None
# Forward solve
    x = np.zeros( n )
    for k in xrange( n - 1 ):
        for i in xrange( k + 1, n ):
            b[p[i]] = b[p[i]] - A[p[i],k] * b[p[k]]
# Backward solve
    for i in xrange( n - 1, -1, -1 ):
        sum = b[p[i]]
        for j in xrange( i + 1, n ):
            sum = sum - A[p[i],j] * x[j]
            x[i] = sum / A[p[i],i]

# All done, return solution vector
    return x

lu(A)
return solve(A,p,v)

def circuit():

A = np.array([[1,0,0,0,0,8,0,0,0],[0,1,0,0,5,0,0,0,0],[0,1,0,0,5,0,0,0,0],[0,0,0,1,-1,1,0,0,0],[0,0,1,0,0,0,1,-1,0],[0,0,1,0,0,0,1,0,-1],[0,1,0,0,-1,0,0,0,1],[1,0,0,0,0,-1,0,1,0],[1,-1,0,1,0,0,0,0,0]])

v = np.array([9,-12,-0.5,0,0,0,0,0,0])
I = solveEqns(A,v)
return I

to solve the 9x9 matrix A at the end. This is one of the easier ones i need to solve so can solve it outside of python to check if the results coming through are accurate.

Im getting a traceback error on line 26 of:

Traceback (most recent call last):

File "<ipython-input-110-6daf773db1e3>", line 1, in <module>
    solveEqns(A,b)

  File "C:/Users/SamMc/Documents/Python Scripts/q6u1510416 v4.py", line 65, in solveEqns
    lu(A)

  File "C:/Users/SamMc/Documents/Python Scripts/q6u1510416 v4.py", line 26, in lu
    r = abs( A[p[i][k]] / s[p[i]] )

TypeError: 'int' object has no attribute '__getitem__'

which i cant figure out why its not pulling through a number from the matrix.

Any help would be greatly appreciated.

Thanks

Sam

Answer 1

you might use gauss elimination via scaled pivoting. the code is shown below.

import numpy as np
def gauss_pivot(a,b,tol=1.0e-12):
    """
    x = gaussPivot(a,b,tol=1.0e-12).
    Solves [a]{x} = {b} by Gauss elimination with
    scaled row pivoting
    """
    a = np.copy(a)
    b = np.copy(b)
    n = len(b)
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix
    # Set up scale factors
    s = np.zeros(n)
    for i in range(n):
        s[i] = max(np.abs(a[i,:])) # find the max of each row
    for k in range(0, n-1): #pivot row
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) # find which row has max item for each col k, and scale by s
        if abs(a[p,k]) < tol: 
            raise Exception("Matrix is singular")
        if p != k: # swap rows if current row does not contain max item with the one contains max item within same col
            a[[k,p+k],:] = a[[p+k, k],:]
            b[k],b[p+k] = b[p+k],b[k] 
            s[k],s[p+k] = s[p+k],s[k]

        # Elimination phase of matrix a
        for i in range(k+1,n):
            if a[i,k] != 0.0: # skip if a(i,k) is already zero
                lam = a [i,k]/a[k,k] 
                a[i,k:n] = a[i,k:n] - lam*a[k,k:n]
                b[i] = b[i] - lam*b[k]

    if abs(a[n-1,n-1]) < tol: 
        raise Exception("Matrix is singular")

    # Back substitution phase, solution is substituted by b
    x = np.zeros_like(b)
    x[n-1] = b[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        x[k] = (b[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k]
    return x

a = np.random.randn(100,100)*10
b = np.random.randn(100)*10

x = gauss_pivot(a,b)

if np.allclose(np.dot(a,x), b) == True:
    print("x is the correct solution")

If you want the code to perform faster you might probably replace x by b, so upon function return b contains the solution. you might also slightly modify elimination phase so elements of matrix a below diagonal are not zeroed, since there are irrelevant during back substitution phase. Therefore, the code becomes as shown below:

import numpy as np
def gauss_pivot(a,b,tol=1.0e-12):
    """
    x = gaussPivot(a,b,tol=1.0e-12).
    Solves [a]{x} = {b} by Gauss elimination with
    scaled row pivoting
    """
    a = np.copy(a)
    b = np.copy(b)
    n = len(b)
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix
    # Set up scale factors
    s = np.zeros(n)
    for i in range(n):
        s[i] = max(np.abs(a[i,:])) # find the max of each row
    for k in range(0, n-1): #pivot row
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) # find which row has max item for each col k, and scale by s
        if abs(a[p,k]) < tol: 
            raise Exception("Matrix is singular")
        if p != k: # swap rows if current row does not contain max item with the one contains max item within same col
            a[[k,p+k],:] = a[[p+k, k],:]
            b[k],b[p+k] = b[p+k],b[k] 
            s[k],s[p+k] = s[p+k],s[k]

        # Elimination phase of matrix a
        for i in range(k+1,n):
            if a[i,k] != 0.0: # skip if a(i,k) is already zero
                lam = a [i,k]/a[k,k] 
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                b[i] = b[i] - lam*b[k]

    if abs(a[n-1,n-1]) < tol: 
        raise Exception("Matrix is singular")

    # Back substitution phase, solution is substituted by b

    b[n-1] = b[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b

To use LU decomposition instead which is more ideal for b containing more than one column, the LU code is shown below

import numpy as np
def lu_decomp(a,tol=1.0e-9):
    a = np.copy(a)
    n = len(a)
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix
    seq = np.arange(n, dtype=int)
    s = np.zeros((n))
    for i in range(n):
        s[i] = max(abs(a[i,:]))
    for k in range(0,n-1):
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) 
        if abs(a[p,k]) < tol: 
            raise Exception("Matrix is singular")
        if p != k:
            a[[k,p+k],:]     = a[[p+k, k],:]
            s[k],s[p+k]      = s[p+k],s[k]
            seq[k], seq[p+k] = seq[p+k],seq[k]

        # Elimination
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                a[i,k] = lam
    return a,seq

def lu_solve(a,b,seq):
    n = len(a)
    x = b.copy()
    for i in range(n):
        x[i] = b[seq[i]]
    # Solution
    for k in range(1,n):
        x[k] = x[k] - np.dot(a[k,0:k],x[0:k])
    x[n-1] = x[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        x[k] = (x[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k]
    return x

a2 = np.random.randn(500,500)*100
b2 = np.random.randn(500,20)*100

a_decomposed, seq = lu_decomp(a2)

x2 = np.zeros_like(b2)

for col in range(b2.shape[1]):
    x2[:,col] = lu_solve(a_decomposed, b2[:, col], seq)


if np.allclose(np.dot(a2,x2), b2) == True:
    print("x2 is the correct solution")

Both methods gives the the output,

Gauss Elimination

x is the correct solution

LU method

x2 is the correct solution

I recommend you use scipy linalg package, from scipy.linalg import solve, lu_factor, lu_solve. They perform way faster for large matrix size. you can use the same code above but annotate them with numba jit so for large matrix the performance is way better.

from numba import jit
@jit
def gauss_pivot(a, b):
    ...
    ...

acknowledgement: codes inspired from the book numerical methods in science and engineering with Python by Prof. Jaan Kiusalaas

https://www.amazon.co.uk/Numerical-Methods-Engineering-Python-3/dp/1107033853/ref=sr_1_1?ie=UTF8&qid=1517845946&sr=8-1&keywords=numerical+method+in+science+and+engineering+with+python