Solve and Simplify

MMA version 12.3 Windows

Has anybody an explanation why "Solve" gives different results if you use "Simplify" on the argument or not.

Here we create 3 equations (it does not matter how they are created):

d = Table[Append[RandomReal[{-1, 1}, {3}], 1], 3];
pl = Cross @@ d;
pt = {x, y, z, 1};
eq1 = pl . pt == 0;
eq2 = ( ((# - pt) . (# - pt)) & /@ d );
eq = {eq1, eq2[[1]] == eq2[[2]], eq2[[2]] == eq2[[3]]};


If I solve this, I get:

sol= Solve[eq, {x, y, z}]


This is not a full solution, x is still free and it does not solve the first equation for arbitrary x. E.g.:

eq /. sol /. x -> 1
(* {False, True, True} *)


However if I simplify the equations I get the full correct solution:

eq = eq // Simplify;
sol=Solve[eq, {x, y, z}][[1]]

eq /. sol
(* {True, True, True} *)

• First, Solve also produces a warning "Solve::svars: Equations may not give solutions for all "solve" variables." Second, Solve may use nonequivalent transformations. Third, use Reduce instead of Solve. Last, but not least: your code is irreproducible because you use RandomReal without Block. – user64494 May 25 at 13:36
• I used "RandomReal" on purpose to show that the error does not depend on a special input. But the point is Solve returns a WRONG result! – Daniel Huber May 25 at 13:40

Well, Solve does not do Simplify on its input automatically. It will not be a good idea for Solve to simplify its input automatically.

When you did simplify, the expression changed just enough to make it give the full answer you saw. (may be due to numerical rounding or cancellation and such, since the input is not exact)

The warning says this basically

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

One way to get same answer as the one you obtained using Simplify but without using Simplify is to convert the input to exact before calling Solve

eqGood = SetPrecision[eq, Infinity];
sol = Solve[eqGood, {x, y, z}] // N


Compare to

sol = Solve[Simplify@eq, {x, y, z}]


Solve is really meant for exact input. Why not just use NSolve? This way, you'll get all the solutions, without simplify and without having to convert input to exact:

 sol = NSolve[eq, {x, y, z}]


• Thank's a lot, your hints are appreciated, but that does not explain why Solve gives a wrong result. – Daniel Huber May 25 at 10:04
SeedRandom[1234];

d = Table[Append[RandomReal[{-1, 1}, {3}], 1], 3];
pl = Cross @@ d;
pt = {x, y, z, 1};
eq1 = pl . pt == 0;
eq2 = (((# - pt) . (# - pt)) & /@ d);
eq = {eq1, eq2[[1]] == eq2[[2]], eq2[[2]] == eq2[[3]]};


Solve is an exact solver and should be given exact input.

sol = Solve[eq // Rationalize[#, 0] &, {x, y, z}] // N

(* {{x -> 0.952153, y -> -0.786618, z -> 0.232846}} *)

eq /. sol[[1]]

(* {True, True, True} *)


Or, for inexact input use NSolve

solN = NSolve[eq, {x, y, z}]

(* {{x -> 0.952153, y -> -0.786618, z -> 0.232846}} *)

eq /. solN[[1]]

(* {True, True, True} *)


Solve without simplification of eq produces an unfinished work, but a correct result. Here are my arguments. Let us execute (BlockRandom guarantees this is reproducible.)

d= BlockRandom[Table[Append[RandomReal[{-1, 1}, {3}], 1], 3]];
pl = Cross @@ d;pt = {x, y, z, 1};eq1 = pl . pt == 0;
eq2 = (((# - pt) . (# - pt)) & /@ d);
eq = {eq1, eq2[[1]] == eq2[[2]], eq2[[2]] == eq2[[3]]}
sol = Solve[eq, {x, y, z}]


{{y -> 0.442888 - 0.456132 x, z -> 0.0787157 - 0.00412639 x}}

and a warning "Solve::svars: Equations may not give solutions for all "solve" variables.". Next,

Simplify[eq /. %[[1]]]


{1. x == 0.602777, True, True}

At last,

Solve[{1. x == 0.602777, True, True}, {x}, Reals]


{{x -> 0.602777}}

Compare

eq


0.324732 - 0.617528 x + 0.281674 y + 0.00254816 z == 0, (0.464599 - x)^2 + (-0.130691 - y)^2 + (-0.399105 - z)^2 == (0.810494 - x)^2 + (0.62591 - y)^2 + (-0.208966 - z)^2, (0.810494 - x)^2 + (0.62591 - y)^2 + (-0.208966 - z)^2 == (0.399552 - x)^2 + (-0.280342 - y)^2 + (0.379474 - z)^2}

with

Simplify[eq]


{1. x == 0.525859 + 0.456132 y + 0.00412639 z, 1. x + 2.18737 y + 0.549701 z == 1.01203, 1. x + 2.2053 y == 0.863986 + 1.43193 z}

Hope this explains the problem.

• @user644943 Thanks for pointing out that the result is unfinished. I do not think that is acceptable for a not too cheap software. that is otherwise great – Daniel Huber May 25 at 15:25
• @DanielHuber: Every command has its limitations. I have demonstrated Solve produces a correct result contrary your claim "Solve returns a WRONG result". – user64494 May 25 at 15:54
• Sorry, I do not agree. The statement: {y -> 0.442888 - 0.456132 x, z -> 0.0787157 - 0.00412639 x} say that the equations are fulfilled for any value of of x. And this is obviously not the case. – Daniel Huber May 25 at 16:47
• @DanielHuber: You don't take into account a warning "Solve::svars: Equations may not give solutions for all "solve" variables." to this result. – user64494 May 25 at 16:50