Solution doesn't solve equation

Question

I'm using Solve[LHS[y]==RHS, y] to generate a solution ruleset, and yet when I apply these rulesets to LHS I'm not getting RHS.

Let me be a little clearer. Sorry that the expressions are so large, but I haven't been able to reproduce this with a simpler set of equations.

justSin is some function of x and y. sinOfArcTan is the same expression but with the substitution y=ArcTan[y2].

justSin = 
 1/2 ((2 Sin[y] Sinh[2 x])  /
      (Cosh[2 x] + Cos[y] Sinh[2 x])  )  ^2

sinOfArctan = 
 1/2 ((2 Sin@ArcTan[y2] Sinh[2 x])  /
      (Cosh[2 x] + Cos@ArcTan[y2] Sinh[2 x])  )   ^2

Now I want to solve the following equation:

rhs = 2 Sinh[2 x]^2 
justSinSol     = FullSimplify@Solve[justSin == rhs, y]
sinOfArctanSol = FullSimplify@Solve[sinOfArctan == rhs, y2]

I get some rules out, the first one's a little complicated but the second is just

{{y2 -> -Csch[2 x]}, {y2 -> Csch[2 x]}}

Now I apply these rules to the original expression:

applyJS = Table[justSin /. justSinSol[[i]], {i, Length[justSinSol]}]
         // FullSimplify // Normal

applySA = Table[sinOfArctan /. sinOfArctanSol[[i]], {i, Length[sinOfArctanSol]}]
        // FullSimplify // Normal

While applyJS gives an answer equivalent to rhs, applySA does not.

(applyJS - rhs) // FullSimplify   (* --> {0,0,0,0}  *)
(applySA - rhs) // FullSimplify   (* --> not zero   *)

I have no idea why this could be the case. I could imagine the introduction of the substitution removing solutions, but I can't see why Mathematica would return a solution which doesn't (seem to) solve the equation.

---

EDIT: As @Feyre points out, the solution is valid, but only for x<0. Why then does Mathematica not return a conditional?

Answer 1

From the docs for Solve

Solve uses non-equivalent transformations to find solutions of transcendental equations and hence it may not find some solutions and may not establish exact conditions on the validity of the solutions found.
...
With Method -> Reduce, Solve uses only equivalent transformations and finds all solutions.

If we follow the guidance given in the docs, we get valid solutions. The condition returned by Solve can be simplified for real x, if desired, as shown below.

sol = FullSimplify@Solve[sinOfArctan == rhs, y2, Method -> Reduce]

realsol = % // FullSimplify[#, 
   TransformationFunctions -> {Automatic, # /. e_Equal :> Reduce[e, x, Reals] &}] &

Solve::useq: The answer found by Solve contains equational condition(s).... A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions.

(*
  {{y2 -> ConditionalExpression[-Csch[2 x], Coth[2 x] + Sqrt[Coth[2 x]^2] == 0]},
   {y2 -> ConditionalExpression[Csch[2 x], Coth[2 x] + Sqrt[Coth[2 x]^2] == 0]},
   {y2 -> ConditionalExpression[-Csch[2 x], Coth[2 x] + Sqrt[Coth[2 x]^2] == 0]},
   {y2 -> ConditionalExpression[Csch[2 x], Coth[2 x] + Sqrt[Coth[2 x]^2] == 0]}}

  {{y2 -> ConditionalExpression[-Csch[2 x], x < 0]},
   {y2 -> ConditionalExpression[Csch[2 x], x < 0]},
   {y2 -> ConditionalExpression[-Csch[2 x], x < 0]},
   {y2 -> ConditionalExpression[Csch[2 x], x < 0]}}
*)

The real solutions are easier to check than the complex ones, because of the condition.

sinOfArctan == rhs /. realsol // Simplify
sinOfArctan == rhs /. realsol // Simplify[#, x < 0] &
(*
  {ConditionalExpression[True, x < 0],  ConditionalExpression[True, x < 0], 
   ConditionalExpression[True, x < 0],  ConditionalExpression[True, x < 0]}

  {True, True, True, True}
*)

SeedRandom[1];
sinOfArctan - rhs /. sol /. x -> # & /@ 
  RandomComplex[{-2 - 2 I, 2 + 2 I}, 5, WorkingPrecision -> 16] // N
(*
  {{Undefined, Undefined, Undefined, Undefined},
   {0.` + 0.` I, 0.` + 0.` I, 0.` + 0.` I, 0.` + 0.` I},
   {Undefined, Undefined, Undefined, Undefined},
   {0.` + 0.` I, 0.` + 0.` I, 0.` + 0.` I, 0.` + 0.` I},
   {Undefined, Undefined, Undefined, Undefined}}
*)

So I guess it's not a bug, but a documented limitation of the default method in Solve.

Answer 2

The solution is valid, for $x<0$

Simplify[(sinOfArctan /. sinOfArctanSol[[2]]) == rhs, 
 Assumptions -> x < 0]

True

Plot[sinOfArctan /. sinOfArctanSol, {x, -Pi, Pi}]
Plot[rhs, {x, -Pi, Pi}]

Note that if we use

sinOfArctanSol = Reduce[sinOfArctan == rhs, y2];

And then run

Simplify[sinOfArctanSol, Assumptions -> {x > 0}]

False

This suggests there is no solution possible for $x\ge0$

I suspect this has to do with trying to fit odd functions into an even function.