Hi Tom, Sorry for late reply, if you are interested in contributing code, then you can send us a Pull Request on Github: https://github.com/sympy/sympy/pulls/
BTW, Wolframalpha also doesn't calculates <http://www.wolframalpha.com/input/?i=LUDecomposition+%28%5B1%2C2%2C3%5D%2C%5B4%2C5%2C6%5D%2C%5B7%2C8%2C9%5D%29> the LUDecomposition of the Matrix you mentioned, since it's singular. Though, your code works for singular matrices, but I am not sure, If the LUDecomposition method in SymPy, should compute LU Decomposition of a singular matrix. Anyways, you can send us a pull Request someone would surely help. *AMiT Kumar* On Thursday, August 6, 2015 at 11:09:06 PM UTC+5:30, Tom H wrote: > > Hi Amit, > > Thanks for the speedy response! If I understand the code > in gauss_jordan_solve correctly, solving A*x=b with multiple right hand > sides requires a call to gauss_jordan_solve for each right hand side b, > which means that the same Gaussian elimination steps are performed on all > but the last column of the augmented matrix [A | b] for each b. Is this > correct? I have attached what is probably the millionth implementation of > Gaussian elimination for computing the LU decomposition of a matrix, in > case someone finds it useful. Let me know if this is not a proper venue for > sharing such code. Thanks again for your input. > Regards, > Tom > > On Tuesday, August 4, 2015 at 10:27:39 PM UTC-7, AMiT Kumar wrote: >> >> Hi Tom, >> >> LUdecomposition() doesn't works for under determined systems, you can use , >> linsolve(), It supports all types of systems: >> >> In [1]: from sympy import * >> >> In [2]: from sympy.solvers.solveset import linsolve >> >> In [3]: A = Matrix([[1,2,3],[4,5,6],[7,8,9]]) >> >> In [4]: b = Matrix([4, 13, 22]) >> >> In [5]: system = (A, b) >> >> In [6]: r, s, t = symbols('r s t') >> >> In [7]: linsolve((A, b), (r, s, t)) >> Out[8]: {(t + 2, -2⋅t + 1, t)} >> >> If you are interested in Matrix form results, you can >> use gauss_jordan_solve() >> >> In [3]: A = Matrix([[1,2,3],[4,5,6],[7,8,9]]) >> >> >> In [4]: b = Matrix([4, 13, 22]) >> >> >> In [5]: sol, params = A.gauss_jordan_solve(b) >> >> In [7]: init_printing() >> >> In [8]: sol >> Out[8]: >> ⎡ τ₀ + 2 ⎤ >> ⎢ ⎥ >> ⎢-2⋅τ₀ + 1⎥ >> ⎢ ⎥ >> ⎣ τ₀ ⎦ >> >> >> >> >> Please note that, both these functions are *not available in any release*, >> you can use them from the *development version*: >> https://github.com/sympy/sympy >> >> -- >> *AMiT Kumar* >> >> On Wednesday, August 5, 2015 at 4:53:15 AM UTC+5:30, Tom H wrote: >>> >>> I came across the following issue when trying to use Sympy to compute an >>> LU decomposition of a matrix. I'd like to determine the number of solutions >>> a system of equations has, for example, if >>> >>> A = sympy.Matrix([[1,2,3],[4,5,6],[7,8,9]]) >>> b = sympy.Matrix([4, 13, 22]) >>> >>> then the system A*x=b has infinitely many solutions, all of which can be >>> written as >>> >>> x0 + t*n >>> >>> where >>> x0 = sympy.Matrix([2,1,0]) >>> n = sympy.Matrix([1,-2,1]) >>> t = sympy.Symbol('t') >>> >>> This solution can be computed from the LU decomposition of A, however >>> the following code fails to compute the factors L and U. >>> >>> >>> import sympy >>> >>> A = sympy.Matrix([[1,2,3],[4,5,6],[7,8,9]]) >>> >>> LU, perm = A.LUdecomposition() >>> Traceback (most recent call last): >>> File "<stdin>", line 1, in <module> >>> File "/Library/Python/2.7/site-packages/sympy/matrices/matrices.py", >>> line 1297, in LUdecomposition >>> combined, p = self.LUdecomposition_Simple(iszerofunc=_iszero) >>> File "/Library/Python/2.7/site-packages/sympy/matrices/matrices.py", >>> line 1341, in LUdecomposition_Simple >>> raise ValueError("No nonzero pivot found; inversion failed.") >>> ValueError: No nonzero pivot found; inversion failed. >>> >>> A is square and noninvertible, which explains the message in the stack >>> trace about the elimination algorithm encountering a zero subcolumn below >>> the pivot position. Regardless, A can still be written as the product of >>> lower triangular L with unit diagonal, and upper triangular U, where >>> L = sympy.Matrix([[1, 0, 0], [4, 1, 0], [7, 2, 1]]) >>> U = sympy.Matrix([[1, 2, 3], [0, -3, -6], [0, 0, 0]]) >>> >>> Is A.LUdecomposition() the wrong function to call if I want to compute >>> the LU decomposition of A in the case where A is noninvertible, or not >>> square? I wrote my own LU decomposition routine to handle these cases, >>> which I'd be happy to contribute, but I'd like to know if this >>> functionality exists in Sympy. >>> >>> Thanks in advance for your input, >>> Tom >>> >> -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/36232ef2-e0bd-42c4-bb9f-ede02922af59%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
