Repository: sympy/sympy
Pull request: linear_coeffs function added to solveset
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests .-->

fixes #15195 
fixes #14987
closes #14999 as a consolidation of work there

#### Brief description of what is fixed or changed

Raise an error when nonlinearity in equations passed to `linear_eq_to_matrix` is encountered.

#### Other comments

This PR arises from the work started by @czgdp1807 to add non-linearity detection into `linear_eq_to_matrix`. The function `linear_coeffs` was written to provide an efficient way to extract linear coefficients without having to simplify an expression: if the expression is linear and has no nonlinear terms appearing in its unsimplified form, the coefficients of the requested symbols will be returned.

```python
    >>> linear_coeffs(x + y*(z*(3*x + 2) + 3), x)
    {0: y*(2*z + 3), x: 3*y*z + 1}
```
It also warns when nonlinear terms are encountered:
```python
    >>> eq = 1/x*(x - 1) + 1/x
    >>> linear_coeffs(eq, x)
    Traceback (most recent call last):
    ...
    ValueError: nonlinear term encountered: 1/x

    >>> linear_coeffs(x*(y + 1) - x*y, x, y)
    Traceback (most recent call last):
    ...
    ValueError: nonlinear term encountered: x*(y + 1)
```

Note: `linsolve` is a higher level function which uses expansion on its arguments before passing them to `linear_eq_to_matrix` and is not (currently) sensitive to removable nonlinearity.

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * `linear_coeffs` added to solveset
  * nonlinear terms encountered by `linear_eq_to_matrix` will raise an error; repeated symbols and non-symbols will also raise an error
  * as a convenience, `linear_eq_to_matrix` now accepts a single Expr as its first argument.
<!-- END RELEASE NOTES -->


Problem statement:
linear_eq_to_matrix proceeds without exception even if equations are nonlinear
Hi,

I would expect the following code to throw an exception (or give me a warning telling me I'm silly), but it does not. `x**2=0` is not a linear equation.

```
from sympy import symbols, linear_eq_to_matrix, linsolve
x = symbols("x")
linear_eq_to_matrix([x**2], [x])
```

Wrong linear_eq_to_matrix if equations not expanded
I noticed that `linear_eq_to_matrix` doesn't yield the correct matrix if the linear equations are not fully expanded. For example
```python
>>> x1, x2, mu = symbols('x1, x2, mu')
>>> equations = Matrix([[mu * (x1 + x2) + mu],
>>>                     [5 * x2 - 1]])
>>> A, b = linear_eq_to_matrix(equations, symbols)
>>> A
```
gives the wrong result `Matrix([[0, 0], [0, 5]])`. This can be fixed by adding `equations = equations.expand()` prior to the `linear_eq_to_matrix` call. Would it make sense to include an `expand` call in the `linear_eq_to_matrix` method itself? Willing to submit the pull request.

Hints from issue discussion:
I would like to work on this issue. Can I proceed working on this?

Sure, anyone can work on any issue.

Some thoughts:

1.  Technically, there is no bug here because the docs say:
> Here `equations` must be a linear system of equations in `symbols`.

2.  But IMHO its still important to have sanity checks in place.  IIRC, there used to be bugs in ODE solvers because of this for example.

3.  It may be important (for other solvers) to have this code as fast as possible.  That might mean no calculating derivatives, for example.

4.  Worst case (if it cannot be done in a fast way), we can have kwarg `checklinearity=True`, and performance-critical library code could set that to False.
I was working on this issue form the past few hours. I have made [PR](https://github.com/sympy/sympy/pull/15000). Please review and kindly give suggestions on improving it further. :-)
> Please review and kindly give suggestions

Well I started too but then you closed it, see my comments in https://github.com/sympy/sympy/pull/14997
@cbm755 I apologize for random actions of mine while making a pull request for solving this issue. Actually, I was unaware of the release notes procedure. Hence, I got confused by the sympy bot. 
I have reopened the [PR](https://github.com/sympy/sympy/pull/14999) which solves the issue to the best of my knowledge. 
I will add the changes you suggested in this PR only. 
Really sorry for the troubles. :-( 
@cbm755 
I have implemented your last suggestion, that is, `check_linearity = False`. This way the standard tests could be passed and if someone wants to confirm the linearity of equation that he or she is passing to the method can do so by setting, `check_linearity = True`.
If it passes all the tests, then I will make appropriate changes in the documentation of the method as well as add appropriate tests. 
Maybe it will be fixed by #14999?  Can you check?  We'd certainly welcome a pull request that adds a test for this, either on top of #14999 or wait until it is merged...
I think you meant
```
A, b = linear_eq_to_matrix(equations, [x1, x2])
```
And it seems this is still an issue on #14999: more specifically, I'm using:
```py
>>> x1, x2, mu = symbols('x1, x2, mu')
>>> equations = [mu * (x1 + x2) + mu, 5 * x2 - 1]
>>> A, b = linear_eq_to_matrix(equations, [x1, x2])
```

#14999 also breaks putting the equations as rows of a matrix: we should probably fix that.


Of course, my bad.
So this will become a ValueError after #14999, certainly better than a wrong answer.

I think your original comment holds: we need to decide if `linear_eq_to_matrix` should automatically expand such things.
#14999 is a convenience function to extract coefficients from a trivially linear system (i.e. fully expanded). If you have a more complicated system of equations, consider either linsolve (which will expand your system and warn if there are persisting non-linear terms) or jacobian:

```python
>>> Matrix([[mu * (x1 + x2) + mu], [5 * x2 - 1]]).jacobian([mu, x1, x2])
Matrix([
[x1 + x2 + 1, mu, mu],
[          0,  0,  5]])
```

Commit message: make sure Eq simplification happens in linsolve

Patch to review:
diff --git a/sympy/physics/continuum_mechanics/beam.py b/sympy/physics/continuum_mechanics/beam.py
index 8d11dbe8594f..387c105d1fda 100644
--- a/sympy/physics/continuum_mechanics/beam.py
+++ b/sympy/physics/continuum_mechanics/beam.py
@@ -1490,13 +1490,13 @@ class Beam3D(Beam):
     >>> b.bc_deflection = [(0, [0, 0, 0]), (l, [0, 0, 0])]
     >>> b.solve_slope_deflection()
     >>> b.slope()
-    [0, 0, l*x*(-l*q + 3*l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I)) + 3*m)/(6*E*I)
-    + q*x**3/(6*E*I) + x**2*(-l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I))
+    [0, 0, l*x*(-l*q + 3*l*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I)) + 3*m)/(6*E*I)
+    + q*x**3/(6*E*I) + x**2*(-l*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I))
     - m)/(2*E*I)]
     >>> b.deflection()
-    [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I))
-    + l*m*x**2/(4*E*I) - l*x**3*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I)
-    + q*x**4/(24*E*I) + l*x*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I)) - q*x**2/(2*A*G), 0]
+    [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I))
+    + l*m*x**2/(4*E*I) - l*x**3*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I)
+    + q*x**4/(24*E*I) + l*x*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I)) - q*x**2/(2*A*G), 0]
 
     References
     ==========
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 17db45e4a7eb..48b8d2a77352 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -472,10 +472,11 @@ def test_Beam3D():
 
     assert b.shear_force() == [0, -q*x, 0]
     assert b.bending_moment() == [0, 0, -m*x + q*x**2/2]
-    assert b.deflection() == [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l**2*q - 2*A*G*l*m
-            + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I)) + l*m*x**2/(4*E*I) - l*x**3*(A*G*l**2*q
-            - 2*A*G*l*m + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I) + q*x**4/(24*E*I)
-            + l*x*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I))
+    assert b.deflection() == [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l*(l*q - 2*m)
+            + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I)) + l*m*x**2/(4*E*I)
+            - l*x**3*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I))
+            - m*x**3/(6*E*I) + q*x**4/(24*E*I)
+            + l*x*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I))
             - q*x**2/(2*A*G), 0]
 
 
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index c0a243b3d92b..3845ce75ab2b 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -19,7 +19,7 @@
 from sympy.core.facts import InconsistentAssumptions
 from sympy.core.numbers import I, Number, Rational, oo
 from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
-                                expand_log)
+                                expand_log, _mexpand)
 from sympy.core.relational import Eq, Ne
 from sympy.core.symbol import Symbol
 from sympy.simplify.simplify import simplify, fraction, trigsimp
@@ -34,7 +34,7 @@
 from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection,
                         Union, ConditionSet, ImageSet, Complement, Contains)
 from sympy.sets.sets import Set
-from sympy.matrices import Matrix
+from sympy.matrices import Matrix, MatrixBase
 from sympy.polys import (roots, Poly, degree, together, PolynomialError,
                          RootOf, factor)
 from sympy.solvers.solvers import (checksol, denoms, unrad,
@@ -42,11 +42,12 @@
 from sympy.solvers.polysys import solve_poly_system
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
-from sympy.utilities.iterables import numbered_symbols
+from sympy.utilities.iterables import numbered_symbols, has_dups
 from sympy.calculus.util import periodicity, continuous_domain
 from sympy.core.compatibility import ordered, default_sort_key, is_sequence
 
 from types import GeneratorType
+from collections import defaultdict
 
 
 def _masked(f, *atoms):
@@ -1806,13 +1807,73 @@ def solvify(f, symbol, domain):
 ###############################################################################
 
 
+def linear_coeffs(eq, *syms):
+    """Return a dictionary with symbols in `syms` keyed to the
+    coefficient of the corresponding symbol in `eq`. The additive
+    constant is returned with key 0.
+
+    Examples
+    ========
+
+    >>> from sympy.solvers.solveset import linear_coeffs
+    >>> from sympy.abc import x, y, z
+
+
+    >>> linear_coeffs(2*x - 1, x, y)
+    {0: -1, x: 2}
+
+    It is not necessary to expand the expression:
+
+    >>> linear_coeffs(x + y*(z*(3*x + 2) + 3), x)
+    {0: y*(2*z + 3), x: 3*y*z + 1}
+
+    But if there are nonlinear or cross terms -- even if they would
+    cancel after simplification-- an error is raised so the situation
+    does not pass silently past the caller's attention:
+
+    >>> eq = 1/x*(x - 1) + 1/x
+    >>> linear_coeffs(eq.expand(), x)
+    {0: 1}
+    >>> linear_coeffs(eq, x)
+    Traceback (most recent call last):
+    ...
+    ValueError: nonlinear term encountered: 1/x
+
+    >>> linear_coeffs(x*(y + 1) - x*y, x, y)
+    Traceback (most recent call last):
+    ...
+    ValueError: nonlinear term encountered: x*(y + 1)
+    """
+    d = defaultdict(list)
+    c, terms = eq.as_coeff_add(*syms)
+    d[0].append(c)
+    for t in terms:
+        m, f = t.as_coeff_mul(*syms)
+        if len(f) != 1:
+            break
+        f = f[0]
+        if f in syms:
+            d[f].append(m)
+        elif f.is_Add:
+            d1 = linear_coeffs(f, *syms)
+            d[0].append(m*d1.pop(0))
+            xf, vf = list(d1.items())[0]
+            d[xf].append(m*vf)
+        else:
+            break
+    else:
+        for k, v in d.items():
+            d[k] = Add(*v)
+        return defaultdict(lambda: S.Zero, d)  # default now zero
+    raise ValueError('nonlinear term encountered: %s' % t)
+
+
 def linear_eq_to_matrix(equations, *symbols):
     r"""
     Converts a given System of Equations into Matrix form.
     Here `equations` must be a linear system of equations in
-    `symbols`. The order of symbols in input `symbols` will
-    determine the order of coefficients in the returned
-    Matrix.
+    `symbols`. Element M[i, j] corresponds to the coefficient
+    of the jth symbol in the ith equation.
 
     The Matrix form corresponds to the augmented matrix form.
     For example:
@@ -1829,84 +1890,104 @@ def linear_eq_to_matrix(equations, *symbols):
      A = [ 3  1  1 ]   b  =   [-6 ]
          [ 2  4  9 ]          [ 2 ]
 
+    The only simplification performed is to convert
+    `Eq(a, b) -> a - b`.
+
+    Raises
+    ======
+
+    ValueError
+        The equations contain a nonlinear term.
+        The symbols are not given or are not unique.
+
     Examples
     ========
 
     >>> from sympy import linear_eq_to_matrix, symbols
-    >>> x, y, z = symbols('x, y, z')
-    >>> eqns = [x + 2*y + 3*z - 1, 3*x + y + z + 6, 2*x + 4*y + 9*z - 2]
-    >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
-    >>> A
-    Matrix([
-    [1, 2, 3],
-    [3, 1, 1],
-    [2, 4, 9]])
-    >>> b
-    Matrix([
-    [ 1],
-    [-6],
-    [ 2]])
-    >>> eqns = [x + z - 1, y + z, x - y]
+    >>> c, x, y, z = symbols('c, x, y, z')
+
+    The coefficients (numerical or symbolic) of the symbols will
+    be returned as matrices:
+
+    >>> eqns = [c*x + z - 1 - c, y + z, x - y]
     >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
     >>> A
     Matrix([
-    [1,  0, 1],
+    [c,  0, 1],
     [0,  1, 1],
     [1, -1, 0]])
     >>> b
     Matrix([
-    [1],
-    [0],
-    [0]])
-
-    * Symbolic coefficients are also supported
-
-    >>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
-    >>> eqns = [a*x + b*y - c, d*x + e*y - f]
-    >>> A, B = linear_eq_to_matrix(eqns, x, y)
-    >>> A
-    Matrix([
-    [a, b],
-    [d, e]])
-    >>> B
-    Matrix([
-    [c],
-    [f]])
-
+    [c + 1],
+    [    0],
+    [    0]])
+
+    This routine does not simplify expressions and will raise an error
+    if nonlinearity is encountered:
+
+    >>> eqns = [
+    ...     (x**2 - 3*x)/(x - 3) - 3,
+    ...     y**2 - 3*y - y*(y - 4) + x - 4]
+    >>> linear_eq_to_matrix(eqns, [x, y])
+    Traceback (most recent call last):
+    ...
+    ValueError:
+    The term (x**2 - 3*x)/(x - 3) is nonlinear in {x, y}
+
+    Simplifying these equations will discard the removable singularity
+    in the first, reveal the linear structure of the second:
+
+    >>> [e.simplify() for e in eqns]
+    [x - 3, x + y - 4]
+
+    Any such simplification needed to eliminate nonlinear terms must
+    be done before calling this routine.
     """
-
     if not symbols:
-        raise ValueError('Symbols must be given, for which coefficients \
-                         are to be found.')
+        raise ValueError(filldedent('''
+            Symbols must be given, for which coefficients
+            are to be found.
+            '''))
 
     if hasattr(symbols[0], '__iter__'):
         symbols = symbols[0]
 
-    M = Matrix([symbols])
-    # initialize Matrix with symbols + 1 columns
-    M = M.col_insert(len(symbols), Matrix([1]))
-    row_no = 1
-
-    for equation in equations:
-        f = sympify(equation)
-        if isinstance(f, Equality):
-            f = f.lhs - f.rhs
-
-        # Extract coeff of symbols
-        coeff_list = []
-        for symbol in symbols:
-            coeff_list.append(f.coeff(symbol))
-
-        # append constant term (term free from symbols)
-        coeff_list.append(-f.as_coeff_add(*symbols)[0])
+    for i in symbols:
+        if not isinstance(i, Symbol):
+            raise ValueError(filldedent('''
+            Expecting a Symbol but got %s
+            ''' % i))
 
-        # insert equations coeff's into rows
-        M = M.row_insert(row_no, Matrix([coeff_list]))
-        row_no += 1
+    if has_dups(symbols):
+        raise ValueError(filldedent('''
+            Symbols must be unique.
+            '''))
+
+    equations = sympify(equations)
+    if isinstance(equations, MatrixBase):
+        equations = list(equations)
+    elif isinstance(equations, Expr):
+        equations = [equations]
+    elif not is_sequence(equations):
+        raise ValueError(filldedent('''
+            Equation(s) must be given as a sequence, Expr,
+            Eq or Matrix.
+            '''))
 
-    # delete the initialized (Ist) trivial row
-    M.row_del(0)
-    A, b = M[:, :-1], M[:, -1:]
+    A, b = [], []
+    for i, f in enumerate(equations):
+        if isinstance(f, Equality):
+            # non-cancelling method of creating expr; we
+            # don't care about the order of the args or
+            # else we would use _unevaluated_Add
+            f = Add._from_args(
+                Add.make_args(f.lhs) +
+                Add.make_args(-f.rhs))
+        co = linear_coeffs(f, *symbols)
+        coeff_list = [co[i] for i in symbols]
+        b.append(-co[0])
+        A.append(coeff_list)
+    A, b = map(Matrix, (A, b))
     return A, b
 
 
@@ -2053,6 +2134,16 @@ def linsolve(system, *symbols):
     >>> linsolve([], x)
     EmptySet()
 
+    * Although removable singularities and nonlinearity
+      will not raise an error, non-linear systems will:
+
+    >>> linsolve([x*(1/x - 1), (y - 1)**2 - y**2 + 1], x, y)
+    {(1, 1)}
+    >>> linsolve([x**2 - 1], x)
+    Traceback (most recent call last):
+    ...
+    ValueError:
+    The term x**2 is nonlinear in {x}
     """
     if not system:
         return S.EmptySet
@@ -2082,19 +2173,9 @@ def linsolve(system, *symbols):
                     symbols for which a solution is being sought must
                     be given as a sequence, too.
                 '''))
-            system = list(system)
-            for i, eq in enumerate(system):
-                try:
-                    # since we are checking it, we might as well take the
-                    # expanded expr that it will give
-                    system[i] = Poly(eq, symbols).as_expr()
-                    if any (degree(eq, sym) > 1 for sym in symbols):
-                        raise PolynomialError
-                except PolynomialError:
-                    raise ValueError(filldedent('''
-                        %s contains non-linear terms in the
-                        variables to be evaluated
-                    ''') % eq)
+            system = [
+                _mexpand(i.lhs - i.rhs if isinstance(i, Eq) else i,
+                recursive=True) for i in system]
             system, symbols, swap = recast_to_symbols(system, symbols)
             A, b = linear_eq_to_matrix(system, symbols)
             syms_needed_msg = 'free symbols in the equations provided'
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 3785645add42..7cf793b987a1 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1019,26 +1019,39 @@ def test_abs_invert_solvify():
 
 def test_linear_eq_to_matrix():
     x, y, z = symbols('x, y, z')
+    a, b, c, d, e, f, g, h, i, j, k, l = symbols('a:l')
+
     eqns1 = [2*x + y - 2*z - 3, x - y - z, x + y + 3*z - 12]
     eqns2 = [Eq(3*x + 2*y - z, 1), Eq(2*x - 2*y + 4*z, -2), -2*x + y - 2*z]
 
-    A, b = linear_eq_to_matrix(eqns1, x, y, z)
+    A, B = linear_eq_to_matrix(eqns1, x, y, z)
     assert A == Matrix([[2, 1, -2], [1, -1, -1], [1, 1, 3]])
-    assert b == Matrix([[3], [0], [12]])
+    assert B == Matrix([[3], [0], [12]])
 
-    A, b = linear_eq_to_matrix(eqns2, x, y, z)
+    A, B = linear_eq_to_matrix(eqns2, x, y, z)
     assert A == Matrix([[3, 2, -1], [2, -2, 4], [-2, 1, -2]])
-    assert b == Matrix([[1], [-2], [0]])
+    assert B == Matrix([[1], [-2], [0]])
 
     # Pure symbolic coefficients
-    from sympy.abc import a, b, c, d, e, f, g, h, i, j, k, l
-    eqns3 = [a*x + b*y + c*z - d, e*x + f*y + g*z - h, i*x + j*y + k*z - l]
+    eqns3 = [a*b*x + b*y + c*z - d, e*x + d*x + f*y + g*z - h, i*x + j*y + k*z - l]
     A, B = linear_eq_to_matrix(eqns3, x, y, z)
-    assert A == Matrix([[a, b, c], [e, f, g], [i, j, k]])
+    assert A == Matrix([[a*b, b, c], [d + e, f, g], [i, j, k]])
     assert B == Matrix([[d], [h], [l]])
 
-    # raise ValueError if no symbols are given
+    # raise ValueError if
+    # 1) no symbols are given
     raises(ValueError, lambda: linear_eq_to_matrix(eqns3))
+    # 2) there are duplicates
+    raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, x, y]))
+    # 3) there are non-symbols
+    raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, 1/a, y]))
+    # 4) a nonlinear term is detected in the original expression
+    raises(ValueError, lambda: linear_eq_to_matrix(Eq(1/x + x, 1/x)))
+
+    assert linear_eq_to_matrix(1, x) == (Matrix([[0]]), Matrix([[-1]]))
+    assert linear_eq_to_matrix(Matrix(
+        [[a*x + b*y - 7], [5*x + 6*y - c]]), x, y) == (
+        Matrix([[a, b], [5, 6]]), Matrix([[7], [c]]))
 
 
 def test_linsolve():
@@ -1133,6 +1146,13 @@ def test_linsolve():
     Eqns = [8*kilo*newton + x + y, 28*kilo*newton*meter + 3*x*meter]
     assert linsolve(Eqns, x, y) == {(-28000*newton/3, 4000*newton/3)}
 
+    # linsolve fully expands expressions, so removable singularities
+    # and other nonlinearity does not raise an error
+    assert linsolve([Eq(x, x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(1/x, 1/x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(y/x, y/x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(x*(x + 1), x**2 + y)], [x, y]) == {(y, y)}
+
 
 def test_solve_decomposition():
     x = Symbol('x')
@@ -1591,6 +1611,27 @@ def test_issue_10397():
     assert solveset(sqrt(x), x, S.Complexes) == FiniteSet(0)
 
 
+def test_issue_14987():
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x**2], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x*(-3/x + 1) + 2*y - a], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x**2 - 3*x)/(x - 3) - 3], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x + 1)**3 - x**3 - 3*x**2 + 7], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x*(1/x + 1) + y], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x + 1)*y], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(1/x, 1/x + y)], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(y/x, y/x + y)], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(x*(x + 1), x**2 + y)], [x, y]))
+
+
 def test_simplification():
     eq = x + (a - b)/(-2*a + 2*b)
     assert solveset(eq, x) == FiniteSet(S.Half)


Code review comments:
diff_hunk: @@ -1829,84 +1890,104 @@ def linear_eq_to_matrix(equations, *symbols):
      A = [ 3  1  1 ]   b  =   [-6 ]
          [ 2  4  9 ]          [ 2 ]
 
+    The only simplification performed is to convert
+    `Eq(a, b) -> a - b`.
+
+    Raises
+    ======
+
+    ValueError
+        The equations contain a nonlinear term.
+        The symbols are not given or are not unique.
+
     Examples
     ========
 
     >>> from sympy import linear_eq_to_matrix, symbols
-    >>> x, y, z = symbols('x, y, z')
-    >>> eqns = [x + 2*y + 3*z - 1, 3*x + y + z + 6, 2*x + 4*y + 9*z - 2]
-    >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
-    >>> A
-    Matrix([
-    [1, 2, 3],
-    [3, 1, 1],
-    [2, 4, 9]])
-    >>> b
-    Matrix([
-    [ 1],
-    [-6],
-    [ 2]])
-    >>> eqns = [x + z - 1, y + z, x - y]
+    >>> c, x, y, z = symbols('c, x, y, z')
+
+    The coefficients (numerical or symbolic) of the symbols will
+    be returned as matrices:
+
+    >>> eqns = [c*x + z - 1 - c, y + z, x - y]
     >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
     >>> A
     Matrix([
-    [1,  0, 1],
+    [c,  0, 1],
     [0,  1, 1],
     [1, -1, 0]])
     >>> b
     Matrix([
-    [1],
-    [0],
-    [0]])
-
-    * Symbolic coefficients are also supported
-
-    >>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
-    >>> eqns = [a*x + b*y - c, d*x + e*y - f]
-    >>> A, B = linear_eq_to_matrix(eqns, x, y)
-    >>> A
-    Matrix([
-    [a, b],
-    [d, e]])
-    >>> B
-    Matrix([
-    [c],
-    [f]])
-
+    [c + 1],
+    [    0],
+    [    0]])
+
+    This routine does not simplify expressions and will raise an error
+    if nonlinearity is encountered:
+
+    >>> eqns = [
+    ...     (x**2 - 3*x)/(x - 3) - 3,
+    ...     y**2 - 3*y - y*(y - 4) + x - 4]
+    >>> linear_eq_to_matrix(eqns, [x, y])
+    Traceback (most recent call last):
+    ...
+    ValueError:
+    The term (x**2 - 3*x)/(x - 3) is nonlinear in {x, y}
+
+    Simplifying these equations will discard the removable singularity
+    in the first, reveal the linear structure of the second:
+
+    >>> [e.simplify() for e in eqns]
+    [x - 3, x + y - 4]
+
+    Any such simplification needed to eliminate nonlinear terms must
+    be done before calling this routine.
     """
-
     if not symbols:
-        raise ValueError('Symbols must be given, for which coefficients \
-                         are to be found.')
+        raise ValueError(filldedent('''
+            Symbols must be given, for which coefficients
+            are to be found.
+            '''))
 
     if hasattr(symbols[0], '__iter__'):
         symbols = symbols[0]
 
-    M = Matrix([symbols])
-    # initialize Matrix with symbols + 1 columns
-    M = M.col_insert(len(symbols), Matrix([1]))
-    row_no = 1
-
-    for equation in equations:
-        f = sympify(equation)
-        if isinstance(f, Equality):
-            f = f.lhs - f.rhs
-
-        # Extract coeff of symbols
-        coeff_list = []
-        for symbol in symbols:
-            coeff_list.append(f.coeff(symbol))
-
-        # append constant term (term free from symbols)
-        coeff_list.append(-f.as_coeff_add(*symbols)[0])
+    for i in symbols:
+        if not isinstance(i, Symbol):
+            raise ValueError(filldedent('''
+            Expecting a Symbol but got %s
+            ''' % i))
 
-        # insert equations coeff's into rows
-        M = M.row_insert(row_no, Matrix([coeff_list]))
-        row_no += 1
+    if has_dups(symbols):
+        raise ValueError(filldedent('''
+            Symbols must be unique.
+            '''))
path: sympy/solvers/solveset.py
text: @user1:
*Please* just use `raise valueError('Symbols must be unique')`.  I don't care how its formatted for users, I don't care if it magically looks better with filldedent.

Ok, maybe this construct is ok for longer error messages (I have not checked) but I refuse to believe its a good idea for short error messages ;-)