# Intro

I have encountered two radical expressions for which I have to find the roots analytically:

  eq1 =  c - e - (2 a x)/(1 + x^2) + Sqrt[2 a - b + d - (2 a)/(1 + x^2)] Sqrt[b - d + (2 a)/(1 + x^2)]


and

eq2 = -((16 a x)/(
1 + x^2)) + (1/((b - d)^2 + (c -
e)^2))(-5 b Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b -
d)^2 + (c - e)^2)] +
5 d Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] + 3 ((b - d)^2 + (c - e)^2) (c - e) +
Sqrt[3 (2 a - b + d) ((b - d)^2 + (c - e)^2) -
5 c Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] +
5 Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c - e)^2)]
e - (16 a ((b - d)^2 + (c - e)^2))/(1 + x^2)]
Sqrt[(-10 a + 3 b - 3 d) ((b - d)^2 + (c - e)^2) +
5 c Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] -
5 Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c - e)^2)]
e + (16 a ((b - d)^2 + (c - e)^2))/(1 + x^2)])


The specific numeric values of my application are:

Data = {a -> 70/1000, b -> 485/10000, c -> 8563/100000, d -> 115/1000,
e -> 148/1000, x -> 0.054170996410719686};


## First equation direct solution

The first one is easy for Mathematica to solve directly using:

Solve[eq1 == 0, x]

Yielding the following result:

{{x -> (-a (-4 c + 4 e) - Sqrt[
a^2 (-4 c + 4 e)^2 -
4 (2 a b + b^2 + c^2 - 2 a d - 2 b d + d^2 - 2 c e +
e^2) (-2 a b + b^2 + c^2 + 2 a d - 2 b d + d^2 - 2 c e +
e^2)])/(2 (-2 a b + b^2 + c^2 + 2 a d - 2 b d + d^2 - 2 c e + e^2))},
{x -> (-a (-4 c + 4 e) + Sqrt[
a^2 (-4 c + 4 e)^2 -
4 (2 a b + b^2 + c^2 - 2 a d - 2 b d + d^2 - 2 c e +
e^2) (-2 a b + b^2 + c^2 + 2 a d - 2 b d + d^2 - 2 c e +
e^2)])/(2 (-2 a b + b^2 + c^2 + 2 a d - 2 b d + d^2 - 2 c e + e^2))}}


Checking the solution with the numerical values of Data gives me the appropriate one.

## Second equation direct solution

Mathematica (12.1) seems not to be able to solve the second one directly:

Solve[eq2 == 0, x]

At least not within reasonable time.

# Solution attempt

It seems clear that the structure of both equations is similar, only being the complexity of the second greater. I am sure that both of them have solutions, it can be checked using the numerical values I have provided.

So I decided to try and manually solve the first one, i.e. apply the required transformations so that it looks like something easily recognizable (polynomial equation), hoping for it to be handy when solving the second one.

## First equation

In order to simplify things, I tried several change of variables, finally finding one suitable for this case:

1/(1 + x^2) -> y, x-> Sqrt[1 - y]/Sqrt[y]


Resulting in:

c - e - 2 a Sqrt[1 - y] Sqrt[y] + Sqrt[2 a - b + d - 2 a y] Sqrt[b - d + 2 a y]


Note: Two equations result depending on positive or negative substitution of x, I will only show the positive one.

Grouping the radicals on one side and squaring:

(c - e)^2 == (2 a Sqrt[1 - y] Sqrt[y] + Sqrt[2 a - b + d - 2 a y] Sqrt[b - d + 2 a y])^2

c^2 - 2 c e + e^2 ==
2 a b - b^2 - 2 a d + 2 b d - d^2 + 8 a^2 y - 4 a b y + 4 a d y -
8 a^2 y^2 +
4 a Sqrt[1 - y] Sqrt[y] Sqrt[2 a - b + d - 2 a y] Sqrt[
b - d + 2 a y]


And grouping and squaring again:

(c^2 - 2 c e + e^2 - 2 a b + b^2 + 2 a d - 2 b d + d^2 - 8 a^2 y +
4 a b y - 4 a d y + 8 a^2 y^2)^2 == (4 a Sqrt[1 - y] Sqrt[y] Sqrt[
2 a - b + d - 2 a y] Sqrt[b - d + 2 a y])^2


Expanding this result and simplifying, the quartic terms are cancelled out and the resulting equation is a quadratic one:

((b - d)^2 + 2 a (-b + d) + (c - e)^2)^2 -
8 a (2 a - b + d) ((b - d)^2 + (c - e)^2) y +
16 a^2 ((b - d)^2 + (c - e)^2) y^2


Which can be solved easily and the result checked against the numerical values.

## Second equation

Due to the complexity of the second equation it is really annoying to do this process manually.

This kind of equations come up quite frequently in the problems I am currently solving, so the question is:

How could I come with a way of automatically solving these equations which all share the same structure but can have different complexity?

Thank you for reading this far.

EDIT 1:

As noted by Bob Hanlon in his answer, both equations are equivalent as they are the solution to almost the same problem through two different paths.

However, using numerical values to solve them is not what I am looking for since the values of these parameters are not constant, they are subject to change in different applications.

Both equations are equivalent due to the topology of this particular problem, but they would not be equivalent in the general case.

EDIT 2:

Is there a way of identifying the radicals on these equations and their coefficients and grouping them together? This would be useful to separate the equation in radical RHS and polynomial LHS and square them independently.

The ideal solution would be to create a function that performs the steps I manually applied to the first equation to any equation with this structure.

Since you know the values of the parameters, substitute the values before solving. This greatly simplifies the problem.

\$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global*"]

data = {a -> 70/1000, b -> 485/10000, c -> 8563/100000, d -> 115/1000,
e -> 148/1000};

eq1 = c - e - (2 a x)/(1 + x^2) +
Sqrt[2 a - b + d - (2 a)/(1 + x^2)] Sqrt[b - d + (2 a)/(1 + x^2)];


The exact solution is

sol1 = Solve[(eq1 /. data) == 0 // Simplify, x]

(* {{x -> (-1782000 + 17 Sqrt[13521852951])/3596381}} *)


Verifying,

eq1 /. data /. sol1[[1]] // FullSimplify

(* 0 *)


The approximate numeric value is

sol1[[1]] // N[#, 17] &

(* {x -> 0.054170996428079429} *)

eq2 = -((16 a x)/(1 +
x^2)) + (1/((b - d)^2 + (c -
e)^2)) (-5 b Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b -
d)^2 + (c - e)^2)] +
5 d Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c - e)^2)] +
3 ((b - d)^2 + (c - e)^2) (c - e) +
Sqrt[3 (2 a - b + d) ((b - d)^2 + (c - e)^2) -
5 c Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] +
5 Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] e - (16 a ((b - d)^2 + (c - e)^2))/(1 +
x^2)] Sqrt[(-10 a + 3 b - 3 d) ((b - d)^2 + (c - e)^2) +
5 c Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] -
5 Sqrt[(-(b - d)^2 - (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] e + (16 a ((b - d)^2 + (c - e)^2))/(1 + x^2)]);


The exact solution is

sol2 = Solve[(eq2 /. data) == 0 // Simplify, x]

(* {{x -> (-1782000 + 17 Sqrt[13521852951])/3596381}} *)


Verifying,

eq2 /. data /. sol2[[1]] // FullSimplify

(* 0 *)


This is identical to the first result

sol1 === sol2

(* True *)


EDIT: The general solution for eq1 can be simplified by restricting the domain to PositiveReals

(sol1a = Solve[eq1 == 0, x, PositiveReals] // FullSimplify)


sol1a[[1]] /. data // Simplify

(* {x -> (-1782000 + 17 Sqrt[13521852951])/3596381} *)

• Thank you for your interest on this question. Indeed, both equations are equivalent as they are the solution to almost the same problem through two different paths. However, using numerical values to solve them is not what I am looking for since the values of these parameters are not constant, they are subject to change in different applications. Commented Jan 27, 2021 at 7:06

A solution for eq2 that works for your data_set is

-((16*a*x)/(1 + x^2)) + (1/((b - d)^2 + (c - e)^2))*(-5*b*
Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)] +

5*d*Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)] +
3*((b - d)^2 + (c - e)^2)*(c - e) +

Sqrt[3*(2*a - b + d)*((b - d)^2 + (c - e)^2) -
5*c*Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)] +

5*Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)]*
e - (16*a*((b - d)^2 + (c - e)^2))/(1 + x^2)]*

Sqrt[(-10*a + 3*b - 3*d)*((b - d)^2 + (c - e)^2) +
5*c*Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)] -

5*Sqrt[(-(b - d)^2 - (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2)]*
e + (16*a*((b - d)^2 + (c - e)^2))/(1 + x^2)]) /.
x -> (Sqrt[(4*
a^2 - (b - d)^2 - (c - e)^2)*((b - d)^2 + (c - e)^2)^5] +
a*((b - d)^2 + (c - e)^2)*(5*b*
Sqrt[-(((b - d)^2 + (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2))] -

5*d*Sqrt[-(((b - d)^2 + (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2))] + 6*b*d*(c - e) -
3*(d^2 + (c - e)^2)*(c - e) + 3*b^2*(-c + e)))/
(((b - d)^2 + (c - e)^2)*(-10*
a^2*((b - d)^2 + (c - e)^2) + ((b - d)^2 + (c - e)^2)^2 +

a*(3*(b - d)*((b - d)^2 + (c - e)^2) +
5*c*Sqrt[-(((b - d)^2 + (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2))] -

5*Sqrt[-(((b - d)^2 + (c - e)^2)*(-4*
a^2 + (b - d)^2 + (c - e)^2))]*e)))


I'm not sure whether this is true for any {a,b,c,d,e}... found it by reducing the coefficient complexity of your equation using intermediate expressions like a1, a2, and so on, then used Solve.

• Thank you for your answer. In this case, I am more interested in creating a way of solving this kind of equations than in the solution to this specific one. Could you create a routine that does that reduction of complexity programmatically? More equations like these come up when solving the system of equations they belong to. Commented Jan 27, 2021 at 11:02
• @GuillermoOliver No. In this case it meant to replace just 3 expressions (not containing the variable x): the first within the large bracket and the other two within the square roots. I don't know why Solve takes so long Commented Jan 27, 2021 at 11:19
• Nice approach. I guess there must be some transformations Solve is trying to apply that significantly increases computing time. Commented Jan 27, 2021 at 18:46

We first define a replacement rule, that disassembles the expressions into structural pieces. Then we write the pure structure using placeholders and solve the bare structure for x. In the results we replace placeholders by the appropriate expressions.

The following replacement pattern will disassemble the large expressions:

repl = (((x1_ : 0) +
x2_ x/(1 + x^2) + (x3_ : 1) ((x4_ : 0) +
Sqrt[x5_ + x6_/(1 + x^2)] Sqrt[x7_ + x8_/(1 + x^2)])) /;
FreeQ[{x1, x2, x3, x4, x5, x6, x7, x8}, x]) :> {x1, x2, x3, x4,
x5, x6, x7, x8}


With the pure structure using placeholders x1,x2,.. we solve the equation for x:

eq = (x1 + x2 x/(1 + x^2) +
x3 (x4 + Sqrt[x5 + x6/(1 + x^2)] Sqrt[x7 + x8/(1 + x^2)])) == 0;
sol= x /. Solve[eq, x];


This gives 4 different lengthy solutions containing placeholders. We now disassemble eq1 or eq2 and replace the placeholders by the actual expressions. E.g. for the first solution and eq1 this would read:

sol[[1]] /. Thread[{x1, x2, x3, x4, x5, x6, x7, x8} -> (eq1 /. repl)]


This gives an expression too large to spell out.

• Thank you very much, this is the kind of approach I was looking for. I was not aware of those functions you are using there and that way of using placeholders. Unfortunately, this results in expressions that are very large, even for the 'simple' equation eq1. Do you think the expressions could be further simplified applying these 'tricks' to the procedure I described in the manual solution? Using placeholders to automate the variable change and radical grouping. Commented Jan 27, 2021 at 18:52
• How to simplify further? Can you reduce the complexity of the structure, so that we have less than 8 parts? Can some symbolic variable be replaced by numbers? An obvious improvement may be to simplify sol, but MMA take so long, that I did not have the patience to wait for it. That's all for the moment, if something more comes to mind, I let you know. Commented Jan 27, 2021 at 20:44
• You may eliminate x8. After the definition of repl add: repl= repl /. x8->(-x6). And in eq replace x8 by -x6 Commented Jan 28, 2021 at 9:12

I have found that the general expression for these equations (at least for these two, will check on more to come) is:

$$x_1+x_2\frac{x}{1+x^2} + \sqrt{x_3+x_2\frac{1}{1+x^2}}\sqrt{x_4-x_2\frac{1}{1+x^2}}$$

Which significantly simplifies the generic expression Daniel Huber proposed.

It can easily be solved:

eq = (x1 + x2 x/(1 + x^2) +
Sqrt[x3 + x2/(1 + x^2)] Sqrt[x4 - x2/(1 + x^2)]) == 0;
sol = x /. Solve[eq, x];


Yielding the generic solution:

$$x = \left\{-\frac{\sqrt{-\left(x_1^2+x_3 (x_2-x_4)\right) \left(x_1^2-x_4 (x_2+x_3)\right)}+x_1 x_2}{x_1^2-x_3 x_4},\frac{\sqrt{-\left(x_1^2+x_3 (x_2-x_4)\right) \left(x_1^2-x_4 (x_2+x_3)\right)}-x_1 x_2}{x_1^2-x_3 x_4}\right\}$$

Using the following code we can calculate the solutions:

replX = (((x1_ + x2_ x/(1 + x^2) +
Sqrt[x3_ + x5_/(1 + x^2)] Sqrt[x4_ + x6_/(1 + x^2)])) /;
FreeQ[{x1, x2, x3, x4, x5, x6}, x]) :> {x1, x2, x3, x4, x5, x6};
sol /. Thread[{x1, x2, x3, x4, x5, x6} -> (eq1 /. replX)]
sol /. Thread[{x1, x2, x3, x4, x5, x6} -> (eq2 /. replX)]


After a FullSimplify, the solutions are, respectively:

{-((-2 a c +
Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c - e)^2)] +
2 a e)/((b - d)^2 + 2 a (-b + d) + (c - e)^2)), (
2 a c + Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] - 2 a e)/((b - d)^2 + 2 a (-b + d) + (c - e)^2)}

{(Sqrt[(4 a^2 - (b - d)^2 - (c - e)^2) ((b - d)^2 + (c - e)^2)^5] +
a ((b - d)^2 + (c -
e)^2) (5 b Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b -
d)^2 + (c - e)^2)] -
5 d Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] + 6 b d (c - e) - 3 (d^2 + (c - e)^2) (c - e) +
3 b^2 (-c + e)))/(((b - d)^2 + (c -
e)^2) (-10 a^2 ((b - d)^2 + (c - e)^2) + ((b - d)^2 + (c -
e)^2)^2 +
a (3 (b - d) ((b - d)^2 + (c - e)^2) +
5 c Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] -
5 Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] e))), (-Sqrt[(4 a^2 - (b - d)^2 - (c -
e)^2) ((b - d)^2 + (c - e)^2)^5] +
a ((b - d)^2 + (c -
e)^2) (5 b Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b -
d)^2 + (c - e)^2)] -
5 d Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] + 6 b d (c - e) - 3 (d^2 + (c - e)^2) (c - e) +
3 b^2 (-c + e)))/(((b - d)^2 + (c -
e)^2) (-10 a^2 ((b - d)^2 + (c - e)^2) + ((b - d)^2 + (c -
e)^2)^2 +
a (3 (b - d) ((b - d)^2 + (c - e)^2) +
5 c Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] -
5 Sqrt[-((b - d)^2 + (c - e)^2) (-4 a^2 + (b - d)^2 + (c -
e)^2)] e)))}


And verify which numeric solution verifies on both equations:

{0.054171, 0.296056}
{-1.04517, 0.054171}


Thank all of you for leading me in the right direction.