-
-
Save mlgill/744d43ba9ec48d6c5afe5a3c3be41bd6 to your computer and use it in GitHub Desktop.
numpy scipy gaussian elimination using LU decomposition with pivoting
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#! /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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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