Skip to content

Instantly share code, notes, and snippets.

@mlgill
Forked from mikofski/Gaussian_elimination.py
Created April 11, 2017 14:07
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mlgill/744d43ba9ec48d6c5afe5a3c3be41bd6 to your computer and use it in GitHub Desktop.
Save mlgill/744d43ba9ec48d6c5afe5a3c3be41bd6 to your computer and use it in GitHub Desktop.
numpy scipy gaussian elimination using LU decomposition with pivoting
#! /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
@sbykl
Copy link

sbykl commented Jan 26, 2021

Hello @mlgill, I am a new Python learner. I am trying to do Gaussian elimination using LU decomposition using Python as well but I am trying to do it with test matrices are stored in the adjacency list (in each row of the file we have three numbers) something like this:
23 3 0.000001370542294
4 4 0.107816040610854
7 4 0.022782277293175
11 4 -0.00921782470662
And file might have 25 or 50 rows.
Can you give me advice on how I could read .txt file and implement the code?
Thanks

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment