# Numerical search of two simultaneous equations for three unknowns

Given the following two equations with three unknowns $$\eta_1, \eta_2$$ and $$\tau$$:

$$\frac{i e^{-i \phi/2}}{\sqrt{2}} = \cos(2 \tau - \eta_1 -\eta_2)\cos(\eta_1-\eta_2) - i\sin(2 \tau - \eta_1-\eta_2)\sin(\eta_1+\eta_2)$$

$$\frac{i e^{i \phi/2}}{\sqrt{2}} = \cos(2 \tau - \eta_1 -\eta_2)\sin(\eta_1-\eta_2) + i\sin(2 \tau - \eta_1-\eta_2)\cos(\eta_1+\eta_2)$$

I'm trying to write code that numerically searches for $$\eta_1, \eta_2$$ and $$\tau$$, each taken from $$[0, 2\pi]$$ but am getting errors (freezes basically). Please advise if there is a simpler way of doing this type of search as compared to what I have. My prototype code is as follows:

(*Define the system of equations for each equation independently*)
equation1[eta1_, eta2_, tau_] :=
I Exp[-I phi/2]/
Sqrt[2] - (Cos[2 tau - eta1 - eta2] Cos[eta1 - eta2] +
I Sin[2 tau - eta1 - eta2] Sin[eta1 + eta2])

equation2[eta1_, eta2_, tau_] :=
I Exp[I phi/2]/Sqrt[2] - (Cos[2 tau - eta1 - eta2] Sin[eta1 - eta2] -
I Sin[2 tau - eta1 - eta2] Cos[eta1 + eta2])

(*Set the value of phi*)
phi = 2.; (*This parameter is fixed*)

(*Range for eta1,eta2,and tau*)
eta1Values = Subdivide[0, 2  Pi, 1000];
eta2Values = Subdivide[0, 2  Pi, 1000];
tauValues = Subdivide[0, 2  Pi, 1000];

(*Iterate over all combinations of eta1,eta2,and tau*)
Do[Do[Do[(*Use NSolve to find the roots for each equation*)
solution1 =
NSolve[{equation1[eta1, eta2, tau] == 0, 0 <= eta1 <= 2 Pi,
0 <= eta2 <= 2 Pi, 0 <= tau <= 2 Pi}, {eta1, eta2, tau},
Reals];
solution2 =
NSolve[{equation2[eta1, eta2, tau] == 0, 0 <= eta1 <= 2 Pi,
0 <= eta2 <= 2 Pi, 0 <= tau <= 2 Pi}, {eta1, eta2, tau},
Reals];
(*Print the solutions*)
If[Length[solution1] > 0,
Print["Solution for equation 1: ", solution1]];
If[Length[solution2] > 0,
Print["Solution for equation 2: ", solution2]];, {tau,
tauValues}], {eta2, eta2Values}], {eta1, eta1Values}]


The solution should exist based on the Solovay-Kitaev decomposition. More specifically based on Eqs. (13)-(15) of the following paper, where I assume $$u:= \frac{i e^{-i \phi/2}}{\sqrt{2}}$$ and $$w:= \frac{i e^{i \phi/2}}{\sqrt{2}}$$.

I think part of the problem is that because I'm are choosing discrete values from [0, 2 Pi, 1000], unlikely we will find perfect roots. Hence the equations will never be satisfied. I might have to add some kind of tolerance like NSolve[{equation1[eta1, eta2, tau] <= 0.001 instead of NSolve[{equation1[eta1, eta2, tau] == 0. I'm not sure if there is anything else dodgy but the code simply freezes when I run it. Any advice on how to improve the code would be most appreciated.

A simple minimization of the square error ought to do it.

phi = 2;

(*Define the system of equations for each equation independently*)
equation1[eta1_, eta2_, tau_] :=
I  Exp[-I  phi/2]/
Sqrt[2] - (Cos[2  tau - eta1 - eta2]  Cos[eta1 - eta2] +
I  Sin[2  tau - eta1 - eta2]  Sin[eta1 + eta2])

equation2[eta1_, eta2_, tau_] :=
I  Exp[I  phi/2]/
Sqrt[2] - (Cos[2  tau - eta1 - eta2]  Sin[eta1 - eta2] -
I  Sin[2  tau - eta1 - eta2]  Cos[eta1 + eta2])

{err, sol} = NMinimize[{
Abs[equation1[eta1, eta2, tau]]^2 +
Abs[equation2[eta1, eta2, tau]]^2,
0 <= eta1 <= 2 Pi, 0 <= eta2 <= 2 Pi, 0 <= tau <= 2 Pi
}, {eta1, eta2, tau}]

{3.42769*10^-18, {eta1 -> 5.49779, eta2 -> 3.14159, tau -> 2.4635}}


Studying those numerical figures, it looks like there is an exact solution for $$\phi=2$$ at $$\eta_1=7 \pi/4 \\ \eta_2=\pi\\ \tau=(4 + 5\pi)/8$$

FullSimplify[equation1[7 Pi/4, Pi, 1/8  (4 + 5  Pi)]] ==
FullSimplify[equation2[7 Pi/4, Pi, 1/8  (4 + 5  Pi)]] == 0

(* True *)


This can be done symbolically in such a way.

phi = 2; TrigReduce[I  Exp[-I  phi/2]/
Sqrt[2] - (Cos[2  tau - eta1 - eta2]  Cos[eta1 - eta2] +
I  Sin[2  tau - eta1 - eta2]  Sin[eta1 + eta2])]


-(1/2) E^-I (-I Sqrt[2] + E^I Cos[2 eta1 - 2 tau] + E^I Cos[2 eta2 - 2 tau] + I Exp[I phi/2]/ Sqrt[2] - (Cos[2 tau - eta1 - eta2] Sin[eta1 - eta2] - I Sin[2 tau - eta1 - eta2] Cos[eta1 + eta2]) I E^ I Cos[2 eta1 + 2 eta2 - 2 tau] - I E^I Cos[2 tau])

TrigReduce[-I*  Exp[I  phi/2]/
Sqrt[2] - (Cos[2  tau - eta1 - eta2]  Sin[eta1 - eta2] -
I  Sin[2  tau - eta1 - eta2]  Cos[eta1 + eta2])]


1/2 (I Sqrt[2] E^I - Sin[2 eta1 - 2 tau] + Sin[2 eta2 - 2 tau] - I Sin[2 eta1 + 2 eta2 - 2 tau] + I Sin[2 tau])

We see the same arguments in the above. In view of it

Reduce[{a == 2  eta1 - 2  tau, b == eta1 - eta2,
c == 2  tau - eta1 - eta2, d == 2  tau}, {eta1, eta2, tau}]


a == b - c && eta1 == b/2 - c/2 + d/2 && eta2 == -(b/2) - c/2 + d/2 && tau == d/2

Now

eq1 = Simplify[ TrigReduce[
I   Exp[-I   phi/2]/
Sqrt[2] - (Cos[2   tau - eta1 - eta2]   Cos[eta1 - eta2] +
I   Sin[2   tau - eta1 - eta2]   Sin[eta1 + eta2])] /. {eta1 ->
b/2 - c/2 + d/2, eta2 -> -(b/2) - c/2 + d/2, tau -> d/2,
a -> b - c}] == 0;
eq2 = Simplify[TrigReduce[
I  Exp[I  phi/2]/
Sqrt[2] - (Cos[2  tau - eta1 - eta2]  Sin[eta1 - eta2] -
I  Sin[2  tau - eta1 - eta2]  Cos[eta1 + eta2])] /. {eta1 ->
b/2 - c/2 + d/2, eta2 -> -(b/2) - c/2 + d/2, tau -> d/2,
a -> b - c}] == 0;
FindInstance[eq1 && eq2 && 0 < b/2 - c/2 + d/2 < 2*Pi &&
0 < -(b/2) - c/2 + d/2 < 2*Pi && 0 < d/2 < 2*Pi, {b, c, d}]


The latest command produces a long output and a warning .

FullSimplify[%]


{{b -> -(\[Pi]/4), c -> 2 (\[Pi] + ArcSin[Sqrt[1/2 (1 - Sin[1])]]), d -> 2 (\[Pi] + ArcTan[(2 (1 + Sqrt[2]) Cos[ 1] (Cos[1] - Sin[1]) (-198 + 140 Sqrt[2] + (-140 + 99 Sqrt[2]) Cos[1] + 140 Sin[1] - 99 Sqrt[2] Sin[1]) Sqrt[-1 + 2/( 1 + Sin[1])])/((-2 + Sqrt[2]) (-4 + 3 Sqrt[2] + (-6 + 4 Sqrt[2]) Cos[1])^2 (-1 + Sin[1]))])}}

FullSimplify[b/2 - c/2 + d/2 /. %]


{-(\[Pi]/8) - ArcSin[Sqrt[1/2 (1 - Sin[1])]] + ArcTan[(-140 + 99 Sqrt[2] + 2 (99 - 70 Sqrt[2]) Sin[1])/(-140 + 99 Sqrt[2] + 2 (-99 + 70 Sqrt[2]) Cos[1])]}

Similarly with eta2 and tau.