Gaussian Elimination using LU decomposition with NumPy and SciPy

This is a fast way to understand LU decomposition and Gaussian elimination

	#! /usr/bin/env python
	"""
	Solve linear system using LU decomposition and Gaussian elimination
	"""
	import numpy as np
	from scipy.linalg import lu, inv
	def gausselim(A,B):
	"""
	Solve Ax = B using Gaussian elimination and LU decomposition.
	A = LU decompose A into lower and upper triangular matrices
	LUx = B substitute into original equation for A
	Let y = Ux and solve:
	Ly = B --> y = (L^-1)B solve for y using "forward" substitution
	Ux = y --> x = (U^-1)y solve for x using "backward" substitution
	:param A: coefficients in Ax = B
	:type A: numpy.ndarray of size (m, n)
	:param B: dependent variable in Ax = B
	:type B: numpy.ndarray of size (m, 1)
	"""
	# LU decomposition with pivot
	pl, u = lu(A, permute_l=True)
	# forward substitution to solve for Ly = B
	y = np.zeros(B.size)
	for m, b in enumerate(B.flatten()):
	y[m] = b
	# skip for loop if m == 0
	if m:
	for n in xrange(m):
	y[m] -= y[n] * pl[m,n]
	y[m] /= pl[m, m]
	# backward substitution to solve for y = Ux
	x = np.zeros(B.size)
	lastidx = B.size - 1 # last index
	for midx in xrange(B.size):
	m = B.size - 1 - midx # backwards index
	x[m] = y[m]
	if midx:
	for nidx in xrange(midx):
	n = B.size - 1 - nidx
	x[m] -= x[n] * u[m,n]
	x[m] /= u[m, m]
	return x
	if __name__ == '__main__':
	x = gausselim(np.array([[3, 2], [1, -4]]), np.array([[5], [10]]))
	print x

view raw Gaussian_elimination.py hosted with ❤ by GitHub