Gauss_elimination

gauss.

## module gaussElimin
''' x = gaussElimin(a,b).
    Solve[a]{b} = {x} by gauss elimination.
'''

import numpy as np 

def gaussElimin(a,b):
  n = len(b)
 #Elmination Phase
  for k in range(0,n-1):
    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]
        b[i] = b[i] - lam*b[k]
 #Back substitution
  for k in range(n-1,-1,-1):
    b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
  return b
  

#!/user/bin/python
## example2_4
  
def vandermode(v):
  n = len(v)
  a = np.zeros((n,n))
  for j in range(n):
    a[:,j] = v**(n-j-1)
  return a
v = np.array([1.0,1.2,1.4,1.6,1.8,2.0])
b = np.array([0.0,1.0,0.0,1.0,0.0,1.0])
a = vandermode(v)
aOrig = a.copy()  # Save original matrix
bOrig = b.copy()  # and the constant vector
x = gaussElimin(a,b)
det = np.prod(np.diagonal(a))
print('x=\n',x)
print('\ndet=',x)
print('\nCheck result: [a]{x}-b =\n',np.dot(aOrig, x)-bOrig)
input("\nPress return to exit")