Repository: sympy/sympy
Pull request: fix(solvers): split separable systems of equations properly
<!-- 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://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #21882

#### Brief description of what is fixed or changed



#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

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
    * A bug when splitting separable systems of equations in solve was fixed.
<!-- END RELEASE NOTES -->


Problem statement:
Incorrect solutions given by `solve` 
The following code tries to find the solutions of a system of (multivariate quadratic) equations. However, the single solution returned by `solve` is not correct.
It seems like `f = 0` should be part of the solution but is not.

```python
from sympy import solve, sympify

equations = [
    "-k*a + b + 5*f/6 + 2*c/9 + 5*d/6 + 4*a/3",
    "-k*f + 4*f/3 + d/2",
    "-k*d + f/6 + d",
    "13*b/18 + 13*c/18 + 13*a/18",
    "-k*c + b/2 + 20*c/9 + a",
    "-k*b + b + c/18 + a/6",
    "5*b/3 + c/3 + a",
    "2*b/3 + 2*c + 4*a/3",
    "-g"
]
equations = [sympify(e) for e in equations]
unknowns = {s for e in equations for s in e.free_symbols}

solutions = solve(equations, unknowns, dict=True)
# Should all be zero but are not:
for solution in solutions:
    for e in equations:
        print(e.subs(solution).simplify())

```

Output:
```
0
f*(-12*k**2 + 28*k - 15)/(12*(k - 1))
0
65*f*(5 - 6*k)/(108*(3*k**2 - 5*k + 2))
5*f*(6*k - 5)*(9*k - 17)/(36*(k - 1)*(3*k - 2))
0
0
35*f*(5 - 6*k)/(18*(3*k**2 - 5*k + 2))
0
```

Expected Output:
```
0
0
0
0
0
0
0
0
0
```

The bug exists in sympy versions `1.8` and `1.7` but NOT in version `1.6`.

In the meantime, is there a way to circumvent this problem until the issue is fixed? We thought about calling `solve` again on the equations of the output but are not sure if that is the best way.

Hints from issue discussion:
I've bisected this to 67fbbe670dbf2e77b4fdebf6493ee4f32f9c8a79 from #18844.

A simple diff that fixes this is:
```diff
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index fd9e6a5fc8..aea9652aa2 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1737,7 +1737,7 @@ def _solve_system(exprs, symbols, **flags):
     if not exprs:
         return []
 
-    if flags.pop('_split', True):
+    if flags.pop('_split', False):
         # Split the system into connected components
         V = exprs
         symsset = set(symbols)
```
However it should be possible to split the system so there must be a bug somehow that was introduced by the changes in that commit.

This is a better fix but I'm not quite sure why it works:
```diff
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index fd9e6a5fc8..dbe97dd12a 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1757,18 +1757,7 @@ def _solve_system(exprs, symbols, **flags):
                 subsyms = set()
                 for e in subexpr:
                     subsyms |= exprsyms[e]
-                subsyms = list(sorted(subsyms, key = lambda x: sym_indices[x]))
-                # use canonical subset to solve these equations
-                # since there may be redundant equations in the set:
-                # take the first equation of several that may have the
-                # same sub-maximal free symbols of interest; the
-                # other equations that weren't used should be checked
-                # to see that they did not fail -- does the solver
-                # take care of that?
-                choices = sift(subexpr, lambda x: tuple(ordered(exprsyms[x])))
-                subexpr = choices.pop(tuple(ordered(subsyms)), [])
-                for k in choices:
-                    subexpr.append(next(ordered(choices[k])))
+                subsyms = sorted(subsyms, key=default_sort_key)
                 flags['_split'] = False  # skip split step
                 subsol = _solve_system(subexpr, subsyms, **flags)
                 if not isinstance(subsol, list):
```
Apparently that comment was committed by me in 67fbbe670dbf2e77b4fdebf6493ee4f32f9c8a79 but the code in that commit doesn't look like something that I would write and I don't understand what the comment means...

@smichr do you understand what this is doing?
> This is a better fix but I'm not quite sure why it works:

That diff gives three test failures in the solvers tests:
```
test_issue_5849
test_issue_5849_matrix
test_issue_18208
```
I think that in the the #5849 case the tests might just be wrong. #18208 seems related and made a change to the ordering here. Actually seeing that a better fix here is:
```diff
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index fd9e6a5fc8..acb12eb421 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1758,17 +1758,6 @@ def _solve_system(exprs, symbols, **flags):
                 for e in subexpr:
                     subsyms |= exprsyms[e]
                 subsyms = list(sorted(subsyms, key = lambda x: sym_indices[x]))
-                # use canonical subset to solve these equations
-                # since there may be redundant equations in the set:
-                # take the first equation of several that may have the
-                # same sub-maximal free symbols of interest; the
-                # other equations that weren't used should be checked
-                # to see that they did not fail -- does the solver
-                # take care of that?
-                choices = sift(subexpr, lambda x: tuple(ordered(exprsyms[x])))
-                subexpr = choices.pop(tuple(ordered(subsyms)), [])
-                for k in choices:
-                    subexpr.append(next(ordered(choices[k])))
                 flags['_split'] = False  # skip split step
                 subsol = _solve_system(subexpr, subsyms, **flags)
                 if not isinstance(subsol, list):
```
I'm still not sure what the removed code does. With that diff though `test_issue_18208` passes but the other two for #5849 still fail. 



Fix plus update to tests:
```diff
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index fd9e6a5fc8..acb12eb421 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1758,17 +1758,6 @@ def _solve_system(exprs, symbols, **flags):
                 for e in subexpr:
                     subsyms |= exprsyms[e]
                 subsyms = list(sorted(subsyms, key = lambda x: sym_indices[x]))
-                # use canonical subset to solve these equations
-                # since there may be redundant equations in the set:
-                # take the first equation of several that may have the
-                # same sub-maximal free symbols of interest; the
-                # other equations that weren't used should be checked
-                # to see that they did not fail -- does the solver
-                # take care of that?
-                choices = sift(subexpr, lambda x: tuple(ordered(exprsyms[x])))
-                subexpr = choices.pop(tuple(ordered(subsyms)), [])
-                for k in choices:
-                    subexpr.append(next(ordered(choices[k])))
                 flags['_split'] = False  # skip split step
                 subsol = _solve_system(subexpr, subsyms, **flags)
                 if not isinstance(subsol, list):
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index ef5a4176c7..9c9ea50e81 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -1316,6 +1316,11 @@ def test_issue_5114_solvers():
 
 
 def test_issue_5849():
+    #
+    # XXX: This system does not have a solution for most values of the
+    # parameters. Generally solve returns the empty set for systems that are
+    # generically inconsistent.
+    #
     I1, I2, I3, I4, I5, I6 = symbols('I1:7')
     dI1, dI4, dQ2, dQ4, Q2, Q4 = symbols('dI1,dI4,dQ2,dQ4,Q2,Q4')
 
@@ -1332,21 +1337,24 @@ def test_issue_5849():
     )
 
     ans = [{
-    I1: I2 + I6,
-    dI1: -4*I2 - 4*I3 - 4*I5 - 10*I6 + 24,
-    I4: -I5 + I6,
-    dQ4: -I5 + I6,
-    Q4: 3*I5/2 - I6/2 - dI4/2,
+    I1: I2 + I3,
+    dI1: -4*I2 - 8*I3 - 4*I5 - 6*I6 + 24,
+    I4: I3 - I5,
+    dQ4: I3 - I5,
+    Q4: -I3/2 + 3*I5/2 - dI4/2,
     dQ2: I2,
     Q2: 2*I3 + 2*I5 + 3*I6}]
 
     v = I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4
     assert solve(e, *v, manual=True, check=False, dict=True) == ans
-    assert solve(e, *v, manual=True) == ans[0]
+    assert solve(e, *v, manual=True, check=False) == ans[0]
+    assert solve(e, *v, manual=True) == []
+    assert solve(e, *v) == []
+
     # the matrix solver (tested below) doesn't like this because it produces
     # a zero row in the matrix. Is this related to issue 4551?
     assert [ei.subs(
-        ans[0]) for ei in e] == [-I3 + I6, I3 - I6, 0, 0, 0, 0, 0, 0, 0]
+        ans[0]) for ei in e] == [0, 0, I3 - I6, -I3 + I6, 0, 0, 0, 0, 0]
 
 
 def test_issue_5849_matrix():
@@ -1373,14 +1381,7 @@ def test_issue_5849_matrix():
         2*I3 + 2*I5 + 3*I6 - Q2,
         I4 - 2*I5 + 2*Q4 + dI4
     )
-    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == {
-    I1: I2 + I6,
-    dI1: -4*I2 - 4*I3 - 4*I5 - 10*I6 + 24,
-    I4: -I5 + I6,
-    dQ4: -I5 + I6,
-    Q4: 3*I5/2 - I6/2 - dI4/2,
-    dQ2: I2,
-    Q2: 2*I3 + 2*I5 + 3*I6}
+    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == []
 
 
 def test_issue_5901():
```
I've opened a PR for this: #21883
@oscarbenjamin We tried out your solution and it seems to resolve the issue mentioned at the beginning but introduces a new problem: Renaming the symbols appearing in the equations changes the solutions (also the number of solutions) we get. For example the following code solves two identical systems then just contains the symbol "_k" instead of "k":

```python
from sympy import solve, sympify, Symbol

equations = [
    "-k*a + b + 5*f/6 + 2*c/9 + 5*d/6 + 4*a/3",
    "-k*f + 4*f/3 + d/2",
    "-k*d + f/6 + d",
    "13*b/18 + 13*c/18 + 13*a/18",
    "-k*c + b/2 + 20*c/9 + a",
    "-k*b + b + c/18 + a/6",
    "5*b/3 + c/3 + a",
    "2*b/3 + 2*c + 4*a/3",
    "-g"
]
equations = [sympify(e) for e in equations]
unknowns = {s for e in equations for s in e.free_symbols}
solutions = solve(equations, unknowns, dict=True)
print(solutions)

equations = [e.subs({Symbol("k"): Symbol("_k")}) for e in equations]
unknowns = {s for e in equations for s in e.free_symbols}
solutions = solve(equations, unknowns, dict=True)
print(solutions)

```

Output:
```
[{a: 0, d: 0, b: 0, f: 0, c: 0, g: 0}, {a: 0, d: -f, b: 0, k: 5/6, c: 0, g: 0}, {a: -2*c, d: 0, b: c, f: 0, k: 13/18, g: 0}, {a: 0, d: 0, b: 0, k: 3/2, c: 0, f: 0, g: 0}]
[{a: 0, _k: 5/6, b: 0, d: -f, c: 0, g: 0}, {a: 0, _k: 3/2, b: 0, d: 0, c: 0, f: 0, g: 0}]
```
> We tried out your solution and it seems to resolve the issue mentioned at the beginning but introduces a new problem:

When you say that this introduces a new problem what do you mean?

I get the same output you showed when testing with master so I don't think it is caused by the change I suggested but is rather a different issue with `solve`. With master:
```
[{a: 0, d: 0, b: 0, f: 0, c: 0, g: 0}, {a: 0, d: -f, b: 0, k: 5/6, c: 0, g: 0}, {a: -2*c, d: 0, b: c, f: 0, k: 13/18, g: 0}, {a: 0, d: 0, b: 0, k: 3/2, c: 0, f: 0, g: 0}]
[{a: 0, _k: 5/6, b: 0, d: -f, c: 0, g: 0}, {a: 0, _k: 3/2, b: 0, d: 0, c: 0, f: 0, g: 0}]
```
With master I get the following output:
```
[{a: 5*f*(6*k - 5)/(12*(k - 1)*(3*k - 2)), b: 0, c: -5*f*(6*k - 5)/(4*(k - 1)*(3*k - 2)), d: f/(6*(k - 1)), g: 0}]
[{_k: (6*d + f)/(6*d), a: 5*d*(d + f)/(2*d + f), b: 0, c: -15*d*(d + f)/(2*d + f), g: 0}]
```
One solution in both cases, although incorrect solutions in both cases for reasons this issue is about, but at least they have the same cardinality.
> With master I get the following output:

Oh, yes I get that now. I must have gotten confused and thought I was testing master when I wasn't.

There's clearly something fishy going on somewhere but it shouldn't be to do with the diff I showed. Those outputs on master are both incorrect. Let's try to work out what the correct answer should be from the Groebner basis:
```python
In [7]: for eq in equations: print(eq)
-a*k + 4*a/3 + b + 2*c/9 + 5*d/6 + 5*f/6
d/2 - f*k + 4*f/3
-d*k + d + f/6
13*a/18 + 13*b/18 + 13*c/18
a + b/2 - c*k + 20*c/9
a/6 - b*k + b + c/18
a + 5*b/3 + c/3
4*a/3 + 2*b/3 + 2*c
-g

In [8]: for eq in groebner(equations).exprs: print(eq)
a + 2*c
b - c
c*f
c*k - 13*c/18
d + f
f*k - 5*f/6
g
```
So we can see that `g=0`, `d=-f`, `b=c` and `a=-2*c`. The remaining unknowns are `c`, `f` and `k` and the remaining equations are
```python
c*f
c*(k - 13/18)
f*(k - 5/6)
```
Now either `c=0` or `f=0`:

Case 1: `c=0`. The final equation is
```python
f*(k - 5/6)
```
and either `f=0` and `k` is anything (case 1a) or `k=5/6` and `f` is anything (case 1b)

Case 2: `f=0`. The final equation is
```python
c*(k - 13/18)
```
and either `c=0` and `k` is anything (case 2a) or `k=13/18` and `c` is anything (case 2b).

The 4 cases given by solve with the diff above are:
```python
{a: 0, d: 0, b: 0, f: 0, c: 0, g: 0}    # case 1a and 2a (same)
{a: 0, d: -f, b: 0, k: 5/6, c: 0, g: 0}    # case 1b
{a: -2*c, d: 0, b: c, f: 0, k: 13/18, g: 0}  # case 2b
{a: 0, d: 0, b: 0, k: 3/2, c: 0, f: 0, g: 0}    # other
```
Now the "other" case is a valid solution but it is in fact a special case of the first solution which has the same zeros but allows `k` to be arbitrary - that case could be eliminated without losing anything from the full solution set since it is implied by the first case.

So the correct result to return should be:
```python
{a: 0, d: 0, b: 0, f: 0, c: 0, g: 0}    # case 1a and 2a (same)
{a: 0, d: -f, b: 0, k: 5/6, c: 0, g: 0}    # case 1b
{a: -2*c, d: 0, b: c, f: 0, k: 13/18, g: 0}  # case 2b
```
Each of these is a valid solution. None are redundant since each solution allows a different symbol ('k', 'f', 'c') to be anything while the other solutions require a particular value for the same symbol. Also I think that the argument above shows that this is a comprehensive set of solutions.

So the diff I showed means that we get correct (although redundant) output for the equations with `k` but not for `k_` so there must be a bug somewhere to do with ordering the symbols. I'm not sure why the output is incorrect on master or why the diff above changes this but there must be a bug.
This seems to make the results for `k` and `k_` consistent:
```diff
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 17540efcc3..a7a6f5beba 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1886,7 +1886,7 @@ def _ok_syms(e, sort=False):
             rv = (e.free_symbols - solved_syms) & legal
             if sort:
                 rv = list(rv)
-                rv.sort(key=default_sort_key)
+                rv.sort(key=lambda s: s.name.replace('_', ''))
             return rv
 
         solved_syms = set(solved_syms)  # set of symbols we have solved for
```
There must be a bug somewhere in this:
https://github.com/sympy/sympy/blob/801b1045d714686fe1baec6ece94815fe43c003a/sympy/solvers/solvers.py#L1890-L1961

That code is called when `solve_poly_system` fails. I'm not sure why `solve_poly_system` would fail here but there's a clearly a bug in the code to handle the failure.

Commit message: fix(solvers): fix doctest after change to solve_system

Patch to review:
diff --git a/doc/src/tutorial/solvers.rst b/doc/src/tutorial/solvers.rst
index f41f1d9f5979..7091ee6a0af1 100644
--- a/doc/src/tutorial/solvers.rst
+++ b/doc/src/tutorial/solvers.rst
@@ -142,9 +142,9 @@ In the ``solveset`` module, the non linear system of equations is solved using
    ``solve`` can be used for such cases:
 
    >>> solve([x**2 - y**2/exp(x)], [x, y], dict=True)
-   ⎡⎧      ⎛-y ⎞⎫  ⎧      ⎛y⎞⎫⎤
-   ⎢⎨x: 2⋅W⎜───⎟⎬, ⎨x: 2⋅W⎜─⎟⎬⎥
-   ⎣⎩      ⎝ 2 ⎠⎭  ⎩      ⎝2⎠⎭⎦
+   ⎡⎧         ____⎫  ⎧        ____⎫⎤
+   ⎢⎨        ╱  x ⎬  ⎨       ╱  x ⎬⎥
+   ⎣⎩y: -x⋅╲╱  ℯ  ⎭, ⎩y: x⋅╲╱  ℯ  ⎭⎦
 
    3. Currently ``nonlinsolve`` is not properly capable of solving the system of equations
    having trigonometric functions.
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index fd9e6a5fc83f..f5c109c0a0f2 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -51,7 +51,7 @@
 from sympy.utilities.lambdify import lambdify
 from sympy.utilities.misc import filldedent
 from sympy.utilities.iterables import (cartes, connected_components,
-    generate_bell, uniq, sift)
+    generate_bell, uniq)
 from sympy.utilities.decorator import conserve_mpmath_dps
 
 from mpmath import findroot
@@ -1758,17 +1758,6 @@ def _solve_system(exprs, symbols, **flags):
                 for e in subexpr:
                     subsyms |= exprsyms[e]
                 subsyms = list(sorted(subsyms, key = lambda x: sym_indices[x]))
-                # use canonical subset to solve these equations
-                # since there may be redundant equations in the set:
-                # take the first equation of several that may have the
-                # same sub-maximal free symbols of interest; the
-                # other equations that weren't used should be checked
-                # to see that they did not fail -- does the solver
-                # take care of that?
-                choices = sift(subexpr, lambda x: tuple(ordered(exprsyms[x])))
-                subexpr = choices.pop(tuple(ordered(subsyms)), [])
-                for k in choices:
-                    subexpr.append(next(ordered(choices[k])))
                 flags['_split'] = False  # skip split step
                 subsol = _solve_system(subexpr, subsyms, **flags)
                 if not isinstance(subsol, list):
@@ -1895,9 +1884,22 @@ def _solve_system(exprs, symbols, **flags):
         # be solved.
         def _ok_syms(e, sort=False):
             rv = (e.free_symbols - solved_syms) & legal
+
+            # Solve first for symbols that have lower degree in the equation.
+            # Ideally we want to solve firstly for symbols that appear linearly
+            # with rational coefficients e.g. if e = x*y + z then we should
+            # solve for z first.
+            def key(sym):
+                ep = e.as_poly(sym)
+                if ep is None:
+                    complexity = (S.Infinity, S.Infinity, S.Infinity)
+                else:
+                    coeff_syms = ep.LC().free_symbols
+                    complexity = (ep.degree(), len(coeff_syms & rv), len(coeff_syms))
+                return complexity + (default_sort_key(sym),)
+
             if sort:
-                rv = list(rv)
-                rv.sort(key=default_sort_key)
+                rv = sorted(rv, key=key)
             return rv
 
         solved_syms = set(solved_syms)  # set of symbols we have solved for
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index ef5a4176c719..9ea118021c9f 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -870,12 +870,12 @@ def test_issue_5132():
         {(log(-sin(2)), -S(2)), (log(sin(2)), S(2))}
     eqs = [exp(x)**2 - sin(y) + z**2, 1/exp(y) - 3]
     assert solve(eqs, set=True) == \
-        ([x, y], {
-        (log(-sqrt(-z**2 - sin(log(3)))), -log(3)),
-        (log(-z**2 - sin(log(3)))/2, -log(3))})
+        ([y, z], {
+        (-log(3), sqrt(-exp(2*x) - sin(log(3)))),
+        (-log(3), -sqrt(-exp(2*x) - sin(log(3))))})
     assert solve(eqs, x, z, set=True) == (
         [x, z],
-        {(log(-z**2 + sin(y))/2, z), (log(-sqrt(-z**2 + sin(y))), z)})
+        {(x, sqrt(-exp(2*x) + sin(y))), (x, -sqrt(-exp(2*x) + sin(y)))})
     assert set(solve(eqs, x, y)) == \
         {
             (log(-sqrt(-z**2 - sin(log(3)))), -log(3)),
@@ -885,12 +885,10 @@ def test_issue_5132():
             (-log(3), -sqrt(-exp(2*x) - sin(log(3)))),
         (-log(3), sqrt(-exp(2*x) - sin(log(3))))}
     eqs = [exp(x)**2 - sin(y) + z, 1/exp(y) - 3]
-    assert solve(eqs, set=True) == ([x, y], {
-        (log(-sqrt(-z - sin(log(3)))), -log(3)),
-            (log(-z - sin(log(3)))/2, -log(3))})
+    assert solve(eqs, set=True) == ([y, z], {
+        (-log(3), -exp(2*x) - sin(log(3)))})
     assert solve(eqs, x, z, set=True) == (
-        [x, z],
-        {(log(-sqrt(-z + sin(y))), z), (log(-z + sin(y))/2, z)})
+        [x, z], {(x, -exp(2*x) + sin(y))})
     assert set(solve(eqs, x, y)) == {
             (log(-sqrt(-z - sin(log(3)))), -log(3)),
             (log(-z - sin(log(3)))/2, -log(3))}
@@ -1316,6 +1314,11 @@ def test_issue_5114_solvers():
 
 
 def test_issue_5849():
+    #
+    # XXX: This system does not have a solution for most values of the
+    # parameters. Generally solve returns the empty set for systems that are
+    # generically inconsistent.
+    #
     I1, I2, I3, I4, I5, I6 = symbols('I1:7')
     dI1, dI4, dQ2, dQ4, Q2, Q4 = symbols('dI1,dI4,dQ2,dQ4,Q2,Q4')
 
@@ -1332,21 +1335,24 @@ def test_issue_5849():
     )
 
     ans = [{
-    I1: I2 + I6,
-    dI1: -4*I2 - 4*I3 - 4*I5 - 10*I6 + 24,
-    I4: -I5 + I6,
-    dQ4: -I5 + I6,
-    Q4: 3*I5/2 - I6/2 - dI4/2,
+    I1: I2 + I3,
+    dI1: -4*I2 - 8*I3 - 4*I5 - 6*I6 + 24,
+    I4: I3 - I5,
+    dQ4: I3 - I5,
+    Q4: -I3/2 + 3*I5/2 - dI4/2,
     dQ2: I2,
     Q2: 2*I3 + 2*I5 + 3*I6}]
 
     v = I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4
     assert solve(e, *v, manual=True, check=False, dict=True) == ans
-    assert solve(e, *v, manual=True) == ans[0]
+    assert solve(e, *v, manual=True, check=False) == ans[0]
+    assert solve(e, *v, manual=True) == []
+    assert solve(e, *v) == []
+
     # the matrix solver (tested below) doesn't like this because it produces
     # a zero row in the matrix. Is this related to issue 4551?
     assert [ei.subs(
-        ans[0]) for ei in e] == [-I3 + I6, I3 - I6, 0, 0, 0, 0, 0, 0, 0]
+        ans[0]) for ei in e] == [0, 0, I3 - I6, -I3 + I6, 0, 0, 0, 0, 0]
 
 
 def test_issue_5849_matrix():
@@ -1373,14 +1379,33 @@ def test_issue_5849_matrix():
         2*I3 + 2*I5 + 3*I6 - Q2,
         I4 - 2*I5 + 2*Q4 + dI4
     )
-    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == {
-    I1: I2 + I6,
-    dI1: -4*I2 - 4*I3 - 4*I5 - 10*I6 + 24,
-    I4: -I5 + I6,
-    dQ4: -I5 + I6,
-    Q4: 3*I5/2 - I6/2 - dI4/2,
-    dQ2: I2,
-    Q2: 2*I3 + 2*I5 + 3*I6}
+    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == []
+
+
+def test_issue_21882():
+
+    a, b, c, d, f, g, k = unknowns = symbols('a, b, c, d, f, g, k')
+
+    equations = [
+        -k*a + b + 5*f/6 + 2*c/9 + 5*d/6 + 4*a/3,
+        -k*f + 4*f/3 + d/2,
+        -k*d + f/6 + d,
+        13*b/18 + 13*c/18 + 13*a/18,
+        -k*c + b/2 + 20*c/9 + a,
+        -k*b + b + c/18 + a/6,
+        5*b/3 + c/3 + a,
+        2*b/3 + 2*c + 4*a/3,
+        -g,
+    ]
+
+    answer = [
+        {a: 0, f: 0, b: 0, d: 0, c: 0, g: 0},
+        {a: 0, f: -d, b: 0, k: S(5)/6, c: 0, g: 0},
+        {a: -2*c, f: 0, b: c, d: 0, k: S(13)/18, g: 0},
+        {a: 0, f: 0, b: 0, k: S(3)/2, c: 0, d: 0, g: 0}
+    ]
+
+    assert solve(equations, unknowns, dict=True) == answer
 
 
 def test_issue_5901():


Code review comments:
diff_hunk: @@ -1373,14 +1379,33 @@ def test_issue_5849_matrix():
         2*I3 + 2*I5 + 3*I6 - Q2,
         I4 - 2*I5 + 2*Q4 + dI4
     )
-    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == {
-    I1: I2 + I6,
-    dI1: -4*I2 - 4*I3 - 4*I5 - 10*I6 + 24,
-    I4: -I5 + I6,
-    dQ4: -I5 + I6,
-    Q4: 3*I5/2 - I6/2 - dI4/2,
-    dQ2: I2,
-    Q2: 2*I3 + 2*I5 + 3*I6}
+    assert solve(e, I1, I4, Q2, Q4, dI1, dI4, dQ2, dQ4) == []
+
+
+def test_issue_21882():
+
+    a, b, c, d, f, g, k = unknowns = symbols('a, b, c, d, f, g, k')
+
+    equations = [
+        -k*a + b + 5*f/6 + 2*c/9 + 5*d/6 + 4*a/3,
+        -k*f + 4*f/3 + d/2,
+        -k*d + f/6 + d,
+        13*b/18 + 13*c/18 + 13*a/18,
+        -k*c + b/2 + 20*c/9 + a,
+        -k*b + b + c/18 + a/6,
+        5*b/3 + c/3 + a,
+        2*b/3 + 2*c + 4*a/3,
+        -g,
+    ]
+
+    answer = [
+        {a: 0, f: 0, b: 0, d: 0, c: 0, g: 0},
+        {a: 0, f: -d, b: 0, k: S(5)/6, c: 0, g: 0},
+        {a: -2*c, f: 0, b: c, d: 0, k: S(13)/18, g: 0},
+        {a: 0, f: 0, b: 0, k: S(3)/2, c: 0, d: 0, g: 0}
path: sympy/solvers/tests/test_solvers.py
text: @user1:
```suggestion
```