# Solve gives duplicate solutions for a particular equation

I stumbled upon the following situation, in which Solve returns the same solution twice (note: the indentation is meant to allow for easy comparison with the variations below)

Solve[(g[0]-f'[x])^2 + (g[0]-f'[x]) (z-z0) y + O[y]^2 == 0, f'[x]]

(* Out: {{f'[x] -> g[0]}, {f'[x] -> g[0]}} *)


When I played around a bit to find the simplest setting in which the duplicate is still given I found that basically any modification removes it:

(* drop square *) Solve[(g[0]-f'[x])   + (g[0]-f'[x]) (z-z0) y + O[y]^2 == 0, f'[x]]
(* drop primes *) Solve[(g[0]-f[x] )^2 + (g[0]-f[x] ) (z-z0) y + O[y]^2 == 0, f[x] ]
(*   z0 -> 0   *) Solve[(g[0]-f'[x])^2 + (g[0]-f'[x])  z     y + O[y]^2 == 0, f'[x]]
(* g[0] -> g0  *) Solve[( g0 -f'[x])^2 + ( g0 -f'[x]) (z-z0) y + O[y]^2 == 0, f'[x]]

(* Out: {{f'[x] -> g[0]}} , etc *)


Indeed, at first sight one might think that this is related to the fact that there's a quadratic equation and the duplicate indicates that the solution has degree two (like perhaps for this question) -- and the duplicate does indeed disappear when the square is removed. Curiously, however, as the above shows, basically any modification does the job.

(Removing the $O(y)^2$ also causes the duplicate answer to disappear, but that really yields a different equation, so I have not included it. Edit: there seems to have been some confusion over this so let me elaborate a bit. I obtained the equation after the Series expansion of $f$ in a second variable, and want to solve it order by order in $y$. In other words, we have two separate equations -- the coefficients at constant and linear order in $y$ -- that I want to solve at the same time. See e.g. power-series solution of differential equations for an example where one typically uses this technique.)

I was just curious whether this is a known bug (I'm using Mathematica 11.0), and whether anyone happens to know its source.

• – Jules Lamers Nov 10 '16 at 15:26
• Even if we usually write $f(x)=\mathcal{O}(x^2)$, it should be $f(x)\in\mathcal{O}(x^2)$. I don't think you equation makes sense. Also I don't know what "$\mathcal{O}(x)^2$" means. – anderstood Nov 10 '16 at 15:27
• @anderstood Could you elaborate on that? (In particular, do you mean in Mathematica or in general?) – Jules Lamers Nov 10 '16 at 15:30
• I'm talking in general; $\mathcal{O}(y^2)$ is not a quantity so I don't understand your equation. See big O notation. – anderstood Nov 10 '16 at 15:33
• Well, "$+O(y^2)$" may be also defined as notation for the usual meaning as given at big O notation, in which case there's no problem there. In any case, that's about semantics and notation in general, which doesn't seem relevant here. If you think your comment is relevant to the behaviour in Mathematica, perhaps you can try writing an answer. – Jules Lamers Nov 10 '16 at 15:42

What you observe is indeed very strange. Small, seemingly irrelevant modifications change whether Solve returns the root once or twice. E.g. just removing the factor z-z0 or changing f'[x] to f[x] causes Solve to return it once instead of twice.

I think I can give at least some insight into this behaviour, and I no longer think that this should be called a bug.

The O[y]^2 at the end of your expression means that it is a SeriesData object, not a usual symbolic expression. It can be converted to a usual expression with Normal.

I was not familiar with the use of SeriesData within Solve. Maybe some others aren't either. So I'll point out this section of the documentation which deals with this:

It tells us that:

• When using series, Solve will equate coefficients one by one and solve the resulting system.

• The system can be constructed using LogicalExpand

Let's see what LogicalExpand gives us so we can see what equations are being solved exactly.

LogicalExpand[(g[0] - f'[x])^2 + (g[0] - f'[x]) (z - z0) y + O[y]^2 == 0]
(*    (z - z0) (g[0] - Derivative[1][f][x]) == 0
&& (g[0] - Derivative[1][f][x])^2 == 0 *)


Removing z-z0:

LogicalExpand[(g[0] - f'[x])^2 + (g[0] - f'[x])  y + O[y]^2 == 0]
(* g[0] - Derivative[1][f][x] == 0 &&
g[0]^2 - 2 g[0] Derivative[1][f][x] + Derivative[1][f][x]^2 == 0 *)


Removing the derivative:

LogicalExpand[(g[0] - f[x])^2 + (g[0] - f[x]) (z - z0) y + O[y]^2 == 0]
(* -z f[x] + z0 f[x] + z g[0] - z0 g[0] == 0 &&
f[x]^2 - 2 f[x] g[0] + g[0]^2 == 0 *)


These systems are all equivalent, but they are stated in different forms.

Essentially they are all equivalent to something of the form:

x-y == 0 && (x-y)^2 == 0


Obviously, x == y is a solution of the system of equations. However, we are solving for a single variable while having a system of two equations, which happen to be consistent with each other.

In this case it doesn't make sense to ask if the root of the system is a single one or a double one. x == y is a single root of x-y == 0 but it is a double root of (x-y)^2 == 0. It depends on which of the two equations we solve.

I imagine that when solving this system, Mathematica drops the first equation, solves the second, then verifies that the solution is consistent with the first one.

But which one should it drop? Very likely, that depends on things such as the name of the variables and their alphabetical order, as well as the precise form of the equations. Essentially it can be considered unpredictable. This phenomenon is fairly common when doing symbolic computations with Mathematica. You will find several posts here on StackExchange which shows that expressions simplify differently if the variables are renamed and thus have a different alphabetical ordering. I think the same kind of thing determines which of the two equations gets dropped and which is solved first.

Much of this is of course guesswork, but it does provide a reasonable explanation of your observations.

Update: Here are two much simpler commands that demonstrate the same phenomenon:

Solve[  (x - 1) == 0 && (x - 1)^2 == 0, x]
(* {{x -> 1}} *)

Solve[a (x - 1) == 0 && (x - 1)^2 == 0, x]
(* {{x -> 1}, {x -> 1}} *)

• Note that Solve[x^2 == 0, x] also returns two solutions {{x -> 0}, {x -> 0}}. – anderstood Nov 10 '16 at 15:54
• @anderstood That's right, but there we only have a quadratic equation, whence two roots, which happen to coincide. As Szabolcs points out this is not the case for the equations here – Jules Lamers Nov 10 '16 at 15:55
• This does indeed remove some unnecessary distraction (viz. the series) from what I found, thanks! (In particular I guess I am happy that Mathematica interprets my input in exactly the way I intended it; I didn't realise that this could be an issue.) But: even if we give the LogicExpanded equations to Solve, doesn't this still seem to be some sort of bug? Solve should "know" there is an equation of degree one, so no multiplicities. – Jules Lamers Nov 10 '16 at 16:00
• ...and thanks for the update! That is certainly a more minimal example of, and gives insight into, what I found. If anything to me that seems to point to this being a bug of Solve, wouldn't you think? – Jules Lamers Nov 10 '16 at 16:23

If you have a differential equation use DSolve rather than Solve

soln = DSolve[(g[0] - f'[x])^2 + (g[0] - f'[x]) (z - z0) y == 0, f, x]

(*  {{f -> Function[{x}, C[1] + x g[0]]},
{f -> Function[{x}, C[1] + x (y z - y z0 + g[0])]}}

f[x] /. soln

(*  {C[1] + x g[0], C[1] + x (y z - y z0 + g[0])}  *)

f'[x] /. soln

(*  {g[0], y z - y z0 + g[0]}  *)

• Yes, but even though my equation involves a derivative, it's really just an algebraic equation for $f'$. Moreover, it's an equation that I want to be solved order by order in $y$; this is why I pointed out that just removing the $O(y)^2$ turns it into another problem. For this new problem both Solve and DSolve give two solutions, which this time are not the identical and there is no issue. – Jules Lamers Nov 10 '16 at 14:52

Edit My answer is wrong: writing something like x+O[x]^2==1 does make sense for MMA. I still leave it because it can help other confused people. What is confusing is that MMA interprets expressions with O[x] as SeriesData, for which the equality == has a special meaning. This may lead to writing equalities which make no sense in maths; for example MMA might solve the equality of two series even if they are not defined in the same neighbourhood.

Your equation (g[0]-f'[x])^2 + (g[0]-f'[x]) (z-z0) y + O[y]^2 == 0 does not make sense because of the big O, which is not a quantity. To make it clear, let's work on a simpler example: x+O[x]^2==1.

Obviously, $\forall \alpha\in\mathbb{R}^*,\ \alpha x^2\in\mathcal{O}(x^2)$. So what should return Solve[x+O[x]^2==1,x]? Of course the solutions of Solve[x+2x^2==1,x] are not the same as those of Solve[x+3x^2==1,x] which are not the same os those of Solve[Normal[x+O[x]^2]==1,x].

Szabolics has given ways of solving coefficient by coefficient (note that LogicalExpand removes the O[y]^2).

• It does make sense for Mathematica. See the link in the middle of my answer. It described what Solve does when presented with such an equation. I was also surprised that this works, and it turns out to have a special meaning within Solve. – Szabolcs Nov 10 '16 at 16:13
• Dear @anderstood, as Szabolics answer shows, my notation certainly makes sense in Mathematica and does exactly what I mean by it. In general the notation $+O(y^2)$ (which Mathematica's O happens to format as O[y]^2) can also be defined to make perfect sense -- see my comment in answer to your comment to my question. (Note that in mathematics, notation such as $x +2\pi \mathbb{Z} := \{x + 2\pi k \mid k \in \mathbb{Z} \}$ or e.g. $x+A$ for $A$ a set is also widely used.) Your example misses the point since there's only one variable, and solving order by order in $x$ does not make sense. – Jules Lamers Nov 10 '16 at 16:16
• @JulesLamers What are the solutions of $1+x+\mathcal{O}(x^2)=0$? – anderstood Nov 10 '16 at 16:22
• @Szabolcs This link says LogicExpand gives the equations obtained by equating corresponding coefficients in the power series. But here we are talking about Solve, which is not used with O in the examples. It says Solve[series1==series2,{a1,a2,…}] _ solve for coefficients in power series_; that's not the form given in Solve[(g[0]-f'[x])^2 + (g[0]-f'[x]) (z-z0) y + O[y]^2 == 0, f'[x]]. – anderstood Nov 10 '16 at 16:24
• @anderstood As an equation that is supposed to hold order by order in $x$ there are none: just like Mathematica treats it, it is really a collection of equations that should be valid at the same time; here $1=0$ appears both at order zero (the coefficient of $x^0$) and one (the coefficient of $x^1$), so the equation is written is not consistent. This is really all related to a rather common technique, see e.g. power series solutions of differential equations – Jules Lamers Nov 10 '16 at 16:25