Repository: sympy/sympy

Issue:
linsolve does not find solutions that solve_linear_system does
For some overdetermined systems, `linsolve` does not find solutions which `solve_linear_system` does. For example:

```python
In [12]: from sympy.solvers.solveset import linsolve
from sympy import *

A, B, J1, J2 = symbols('A B J1 J2')

augmatrix = Matrix([
    [ 2*I*J1,  2*I*J2,-2/J1],
    [-2*I*J2, -2*I*J1, 2/J2],
    [ 0,  2, 2*I/(J1*J2)],
    [2, 0,  0],
])

linsolve(augmatrix, A, B)
solve_linear_system(augmatrix, A, B)

Out[16]: EmptySet()

Out[17]: {A: 0, B: I/(J1*J2)}
```
On the other hand:

```python
In [23]: augmatrix = Matrix([
    [2*I*J2,-2/J1],
    [-2*I*J1, 2/J2],
    [ 2, 2*I/(J1*J2)],
])

linsolve(augmatrix, A)
solve_linear_system(augmatrix, A)

Out[24]: {(I/(J1*J2),)}

Out[25]: {A: I/(J1*J2)}
```

works fine.

I understand that this issue is problematic because `solve_linear_system` is being deprecated in favour of  `linsolve`: #10000.


Discussion:

Patch:
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 96d2a4c295..0a61971a3e 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -4325,7 +4325,7 @@ def gauss_jordan_solve(self, b, freevar=False):
         row, col = aug[:, :-1].shape
 
         # solve by reduced row echelon form
-        A, pivots = aug.rref()
+        A, pivots = aug.rref(simplify=True)
         A, v = A[:, :-1], A[:, -1]
         pivots = list(filter(lambda p: p < col, pivots))
         rank = len(pivots)
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 132c550e37..ae65843752 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -12,8 +12,8 @@
 from sympy.core.numbers import I, Number, Rational, oo
 from sympy.core.function import (Lambda, expand, expand_complex)
 from sympy.core.relational import Eq
+from sympy.simplify.simplify import simplify, fraction, trigsimp
 from sympy.core.symbol import Symbol
-from sympy.simplify.simplify import fraction, trigsimp
 from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
                              acos, asin, atan, acsc, asec, arg,
                              Piecewise, piecewise_fold)
@@ -1210,7 +1210,7 @@ def linsolve(system, *symbols):
     >>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
     >>> eqns = [a*x + b*y - c, d*x + e*y - f]
     >>> linsolve(eqns, x, y)
-    {(-b*(f - c*d/a)/(a*(e - b*d/a)) + c/a, (f - c*d/a)/(e - b*d/a))}
+    {((-b*f + c*e)/(a*e - b*d), (a*f - c*d)/(a*e - b*d))}
 
     * A degenerate system returns solution as set of given
       symbols.
@@ -1275,11 +1275,10 @@ def linsolve(system, *symbols):
         for s in sol:
             for k, v in enumerate(params):
                 s = s.xreplace({v: symbols[free_syms[k]]})
-            solution.append(s)
-
+            solution.append(simplify(s))
     else:
         for s in sol:
-            solution.append(s)
+            solution.append(simplify(s))
 
     # Return solutions
     solution = FiniteSet(tuple(solution))


Test patch:
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index ba77f552fc..7e2b8c389d 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -979,8 +979,7 @@ def test_linsolve():
     A = Matrix([[a, b], [c, d]])
     B = Matrix([[e], [f]])
     system2 = (A, B)
-    sol = FiniteSet((-b*(f - c*e/a)/(a*(d - b*c/a)) + e/a,
-                    (f - c*e/a)/(d - b*c/a)))
+    sol = FiniteSet(((-b*f + d*e)/(a*d - b*c), (a*f - c*e)/(a*d - b*c)))
     assert linsolve(system2, [x, y]) == sol
 
     # Test for Dummy Symbols issue #9667
@@ -996,6 +995,17 @@ def test_linsolve():
     b = Matrix([0, 0, 1])
     assert linsolve((A, b), (x, y, z)) == EmptySet()
 
+    # Issue #10056
+    A, B, J1, J2 = symbols('A B J1 J2')
+    Augmatrix = Matrix([
+        [2*I*J1, 2*I*J2, -2/J1],
+        [-2*I*J2, -2*I*J1, 2/J2],
+        [0, 2, 2*I/(J1*J2)],
+        [2, 0,  0],
+        ])
+
+    assert linsolve(Augmatrix, A, B) == FiniteSet((0, I/(J1*J2)))
+
     # Issue #10121 - Assignment of free variables
     a, b, c, d, e = symbols('a, b, c, d, e')
     Augmatrix = Matrix([[0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0]])