import numpy as np
This function performs LU decomposition of given matrix A.
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 $$
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.
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.
b = np.random.rand(n)
x = LUSolve(L,U,b)
print A.dot(x) - b
[ 0.00000000e+00 0.00000000e+00 2.22044605e-16]