# Solving a two variable two equation trig expression

I am trying to use Solve to get what two angles (A and B) are in a problem I have.

Solve[1/2 (-1 + Cos[2 B] + 2 Sqrt Sqrt[Cos[2 A] + Cos[2 B]] Sin[A]) == x &&
(Sqrt[Cos[2 A] + Cos[2 B]] + Sqrt Sin[A]) Sin[B] == z,
{A, B}]


but when I tun it through Mathematica it just returns an empty set. If it matters, I care about the region where $-\pi/4 < \mathtt{A,B} < \pi/4$ but even adding that constraint in didn't help me get a solution from Mathematica.

• I wanted A and B solved in terms of x and z. Like A -> ArcCos(x/z) Am I doing something wrong in expecting that? – Michael Jun 4 '15 at 20:29
• I don't believe Mma will be able to do that in the general case – Dr. belisarius Jun 4 '15 at 20:33
• One trick you can do is to apply TrigExpand[], make the Weierstrass substitution {Cos[A] :> (1 - u^2)/(1 + u^2), Sin[A] :> 2 u/(1 + u^2)} (and similarly for B), feed the algebraic system to Solve[] (or maybe even GroebnerBasis[]) and then undo the substitution. – J. M.'s technical difficulties Jun 4 '15 at 20:35
• @bel, gotta say, I'm not quite sure why it takes so long. But, what I suggested is a standard tool for dealing with algebraic combinations of trigonometric functions. Maybe Daniel can offer his insight, should he come across this. – J. M.'s technical difficulties Jun 5 '15 at 10:33
• @Guesswhoitis. I tried to understand the reason also, but I hadn't enough spare time – Dr. belisarius Jun 5 '15 at 12:32

It turns out that the trivial substitution {Cos@u-> x, Sin@u ->Sqrt[1-x^2]} is easier to work with:

eqs = {1/2 (-1 + Cos[2 B] + 2 Sqrt Sqrt[Cos[2 A] + Cos[2 B]] Sin[A]),
(Sqrt[ Cos[2 A] + Cos[2 B]] + Sqrt Sin[A]) Sin[B]};
r = MapAll[TrigExpand, #] & /@ eqs;
subs[A_, u_] := {Cos[A] :> u, Sin[A] :> Sqrt[1 - u u]}
r1 = r /. subs[A, u] /. subs[B, v];
sols = Solve[Thread[r1 == {x, z}], {u, v}] // FullSimplify;
Partition[
Plot3D[{ArcCos@u, ArcCos@v} /. #, {x, -1, 1}, {z, -1, 1},
Evaluated -> True, PlotStyle -> {Red, Blue},
PlotRange -> {{-1, 1}, {-1, 1}, Pi/4 {-1, 1}},
ClippingStyle -> None] & /@ sols, 4] // GraphicsGrid You can see that from the eight solutions, only the {4, 8} give results within the desired domain. Moreover, you can check that the result for B is the same in both cases:

sols[[4, 2]] == sols[[8, 2]]
(* True *)
sols[[4, 2]]
(*v -> Sqrt[2 + 2 x - z^2]/Sqrt[2 + 2 x]*)


So the only solution for B is

$$\text{B}(x,z) = \cos ^{-1}\left(\frac{\sqrt{2 x-z^2+2}}{\sqrt{2 x+2}}\right)$$

B[x_,z_]: = ArcCos[Sqrt[2 + 2 x - z^2]/Sqrt[2 + 2 x]]

Plot3D[B[x, z], {x, -1, 1}, {z, -1, 1},
PlotRange -> {{-1, 1}, {-1, 1}, Pi/4 {-1, 1}}, ClippingStyle -> None,
PlotStyle -> Blue] For A we have two possible solutions

A1[x_, z_] := Evaluate[ArcCos@u /. sols[[4, 1]]]
A2[x_, z_] := Evaluate[ArcCos@u /. sols[[8, 1]]] Plot3D[{A1[x, z], A2[x, z]}, {x, -1, 1}, {z, -1, 1},
PlotRange -> {{-1, 1}, {-1, 1}, Pi/4 {-1, 1}}, ClippingStyle -> None,
PlotStyle -> {Red, Yellow}] $$\text{A1}(x,z) = \cos ^{-1}\left(\frac{1}{2} \sqrt{\frac{(x+1) z^2-2 \sqrt{-(x+1)^4 \left(x^2+z^2-1\right)}}{(x+1)^2}+2}\right)$$

$$\text{A2}(x,z) = \cos ^{-1}\left(\frac{1}{2} \sqrt{\frac{(x+1) z^2+2 \sqrt{-(x+1)^4 \left(x^2+z^2-1\right)}}{(x+1)^2}+2}\right)$$

This answer is based on @Guesswhoitis's hints on using the Weierstrass substitution and GroebnerBasis in the comments below the question.
Your problem looks deceivingly easy, but perhaps it is not so. Let's see:

subs[A_, u_] := {Cos[A] :> (1 - u^2)/(1 + u^2), Sin[A] :> 2 u/(1 + u^2)}

eq1 = 1/2 (-1 + Cos[2 B] + 2 Sqrt Sqrt[Cos[2 A] + Cos[2 B]] Sin[A]) == x;
eq2 = (Sqrt[Cos[2 A] + Cos[2 B]] + Sqrt Sin[A]) Sin[B] == z;

r = MapAll[TrigExpand, #] /. subs[A, u] /. subs[B, v] & /@ {eq1, eq2};
s = Solve[Thread[Equal[GroebnerBasis[r, {x, z}, {u, v}], 0]], {u, v}];


Warning note: The expressions in s are pretty weird monsters that take quite a CPU spend to simplify, calculate, etc. Don't try to work with them, just look:

{Depth[s], LeafCount[s], Length@s}
(*{25, 1984657, 16}*)


So we are going to use a compiled version (please note that there are 16 solutions):

fun = Compile[{{x, _Real}, {z, _Real}}, #, "RuntimeOptions" -> "Speed"] & /@
Partition[({u, v} /. s // Flatten), 2];
(* and another function to return a rule {u,v}->{something} *)
rfun = MapThread[Rule, {{u, v}, fun[[#1]][#2, #3]}] &;


In fact we made 16 functions instead of evaluating all solutions in one function. That is because at any point in {{-1,1},{-1,1}} there are eight functions that return Complex solutions and then compiled evaluation will be discarded always for all of them. By compiling each solution separately we get that behavior too, but in a restricted domain for each one.

Now we build a test for the (numerical) validity of each of the 16 solutions:

ClearAll[testSol];
testSol[n_?NumericQ, x0_?NumericQ, z0_?NumericQ] :=
If[Abs[Subtract @@ eq1] < 10^-5 &&  Abs[Subtract @@ eq2] < 10^-5 /.
{A -> 2 ArcTan@u, B -> 2 ArcTan@v} /. rfun[n, x0, z0] /. {x -> x0, z -> z0}, 1, 0,0]


Let's take a look at the domain of each solution:

Partition[
Quiet@ContourPlot[testSol[#,x00,z00], {x00, -1, 1}, {z00, -1, 1}, PlotPoints-> 30] &/@
Range@16, 4] // GraphicsGrid So,"roughly", the two columns on the right and the two on the left give complementary domains for the solutions, and you'll have eight solutions for each point in {{-1,1},{-1,1}}.

We can see that the solutions for A on the following subsets of solutions {{1, 2}, {5, 6}, {9, 10}, {13, 14}}} are equal (on each pair), but the solutions for B aren't:

qq = {Quiet@
Plot3D[{2 ArcTan@u /. rfun[#[], x0, z0], 2  ArcTan@u /. rfun[#[], x0, z0]},
{x0, -1, 1}, {z0, -1, 1}, PlotStyle -> {Red, Blue},
RegionFunction -> (testSol[1, #1, #2] == 1 &),
Evaluated -> True] & /@ {{1, 2}, {5, 6}, {9, 10}, {13, 14}}};

GraphicsGrid@qq But when we perform the same for B, you can see that the solutions are "out of phase". So here we have effectively eight different solutions. Please note that due to the discontinuous nature of the solutions for B the plots take a lot of time to run. Lets see what happens with the other eight {{3, 4}, {7, 8}, {11, 12}, {15, 16}}}. First for A

qq = {Quiet@
Plot3D[{2 ArcTan@u /. rfun[#[], x0, z0], 2  ArcTan@u /. rfun[#[], x0, z0]},
{x0, -1, 1}, {z0, -1, 1}, PlotStyle -> {Red, Blue},
RegionFunction -> (testSol[3, #1, #2] == 1 &),
Evaluated -> True] & /@ {{3, 4}, {7, 8}, {11, 12}, {15, 16}}};

GraphicsGrid@qq We can see that the same behavior is reproduced on this side. The solutions for A are equal in each pair. We then explore the behavior for B with PlotPoints -> 10, MaxRecursion -> 0 because the functions are discontinuous and expensive to evaluate. Anyway we are only interested in classifying them by now. So now we know that our 16 solutions are different.

Furthermore you said that you're only interested in solutions for A and B satisfying $-\pi/4 < \mathtt{A,B} < \pi/4$. Let's see what that means:

ClearAll[testSol1];
testSol1[n_?NumericQ, x0_?NumericQ, z0_?NumericQ] :=
If[Abs[Subtract @@ eq1] < 10^-5 &&   Abs[Subtract @@ eq2] < 10^-5 &&
-Pi/4 <= A <= Pi/4 && -Pi/4 <= B <= Pi/4 /. {A -> 2 ArcTan@u, B -> 2 ArcTan@v} /.
rfun[n, x0, z0] /. {x -> x0, z -> z0}, 1, 0, 0]

Partition[
Quiet@ContourPlot[testSol1[#, x00, z00], {x00, -1, 1}, {z00, -1, 1},
PlotPoints -> 30] & /@ Range, 4] // GraphicsGrid You can see from the plots above that only solutions 10 and 11 seem to satisfy the bounds for A and B. However, it is possible that the rest of the solutions have much smaller compliance domains hindered here by numerical instabilities. I won't go further analyzing that (Post Scriptum: I should have.See my other answer). Let's see only how sols 10 and 11 for A and B look like in their corresponding domains:

p[n_] := Quiet@ Plot3D[{2 ArcTan @u, 2 ArcTan@v} /. rfun[n, x0, z0],
{x0, -1,  1}, {z0, -1, 1},
PlotStyle -> {Red, Blue},
RegionFunction -> (testSol[n, #1, #2] == 1 &),
Evaluated -> True]

{p /@ {10, 11}} // GraphicsGrid so they seem to define a "single patched" continuous solution:

Show @@ Cases[%, Inset[x_, ___] :> x, Infinity] • The last plot looks tantalizing. I wonder if there might be a way to parametrize the set of solutions… – J. M.'s technical difficulties Jun 9 '15 at 22:34
• @Guesswhoitis. It seems that the trivial one subs[A_, u_] := {Cos[A] :> u, Sin[A] :> Sqrt[1 - u u]} works much easier. Going to post it in a while – Dr. belisarius Jun 9 '15 at 23:22
• Good to hear that works here; most of the time, I end up with hanging square roots with that substitution, but maybe I just happen to frequently deal with sets where that fails. – J. M.'s technical difficulties Jun 9 '15 at 23:26
• @Guesswhoitis. I posted another answer :( – Dr. belisarius Jun 10 '15 at 0:32