Gaussian Elimination

In [77]:
import numpy as np

This function performs LU decomposition of given matrix A.

In [86]:
def LU(A):
    n = A.shape[0]
    L = np.zeros((n,n))
    U = np.array(A)
    for i in range(n-1):
        L[i,i] = 1.0
        L[i+1:n,i] = U[i+1:n,i]/U[i,i]
        for j in range(i+1,n):
            U[j,i+1:n] = U[j,i+1:n] - L[j,i] * U[i,i+1:n]
        U[i+1:n,i] = 0.0
    L[n-1,n-1] = 1.0
    return L,U

This performs performs solution of $$ LUx=b $$ in two steps. In first step, solve $$ Ly = b $$ In second step, solve $$ Ux = y $$

In [86]:
def LUSolve(L,U,b):
    n = L.shape[0]
    # solve Ly = b
    y = 0*b
    for i in range(n):
        y[i] = (b[i] - L[i,0:i].dot(y[0:i]))/L[i,i]
    # solve Ux = y
    x = 0*b
    for i in range(n-1,-1,-1):
        x[i] = (y[i] - U[i,i+1:n].dot(x[i+1:n]))/U[i,i]
    return x

Now we test the above function for LU decomposition.

In [87]:
n = 3
A = np.random.rand(n,n)
L,U = LU(A)
print A
print L
print U
print A - L.dot(U)
[[ 0.23000755  0.32161034  0.37253327]
 [ 0.87099231  0.62303385  0.06389315]
 [ 0.24376526  0.92755117  0.03188404]]
[[ 1.          0.          0.        ]
 [ 3.78679872  1.          0.        ]
 [ 1.05981414 -0.98632268  1.        ]]
[[ 0.23000755  0.32161034  0.37253327]
 [ 0.         -0.59483979 -1.34681537]
 [ 0.          0.         -1.69132653]]
[[ 0.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]]

Solve the linear system.

In [88]:
b = np.random.rand(n)
x = LUSolve(L,U,b)
print A.dot(x) - b
[  0.00000000e+00   0.00000000e+00   2.22044605e-16]