# Problems with Solve when constraining solutions to a range

I have the following code

{e1, e2, e3} =
{Sin[rot + 2 Pi (# – 1)/3], -Cos[rot + 2 Pi (# – 1)/3]} & /@ {1, 2, 3};
{c1, c2, c3} =
1 – 3 (Sin[theta] Cos[phi – 2 Pi (# – 1)/3 – rot])^2 & /@ {1, 2,3};
fq = Total[#[[1]] Exp[I q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3, e3}}];
{theta, phi} = {Pi/2, Pi/2};
rot = 0;
q = {qx, qy};
fqReal = Total[#[[1]] Cos[ q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3,e3}}]
fqImag = Total[#[[1]] Sin[ q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3,e3}}]
{qx, qy} /.
Solve[fqReal == 0 && fqImag == 0(*&&-Pi<=qx<=Pi&&-Pi<=qy<=Pi*), {qx, qy}] When I leave the commented-out condition within Solve it works out the solution and gives me the answer in the form of a list of ConditionalExpressions. However when I ask it to also include that condition -Pi<=qx<=Pi&&-Pi<=qy<=Pi I get the error warning "Solve::nsmet: This system cannot be solved with the methods available to Solve." From the physics, I feel quite confident that there is mathematically a solution within this range, so why am I getting the error warning and how do I fix this problem? Perhaps alternatively, a solution would be to take the ConditionalExpressions and somehow filter them using the condition -Pi<=qx<=Pi&&-Pi<=qy<=Pi, however I do not know how to do this. =================      Not every equation can be solved with Solve. In this case: try NSolve, to solve the eqn. numerically – mgamer Jan 4 '15 at 13:24      That's what I figured but unfortunately that doesn't work either 🙁 – Tom Jan 4 '15 at 14:25 ================= 2 Answers 2 ================= The problem with arbitary rot, mentioned in a comment by the OP, is that the equations contain terms like (1 - 3 Cos[Ï€/6 + rot]^2) Cos[qx Cos[Ï€/6 + rot] + qy Sin[Ï€/6 + rot]] That rot appears inside Cos simply as well as inside Sin and Cos which are themselves inside a Cos, and that there are several such terms, makes this a difficult equation to solve symbolically. But we can differentiate the equations with respect to rot and use the solutions when rot = 0 as initial conditions. The equations can them be integrated to give numerical approximations to the solutions. I didn't know what range rot might have so I picked {rot, -Ï€, Ï€} initially. NDSolve emits a message indicating that it should have been Â±0.178. All of the IVPs have this same stopping point, which must be more than a coincidence and might be deducible from the ODE (if someone is curious). Note: Solve can call Reduce, as in Dr. Wolfgang Hintze's answer, with the option Method -> Reduce.

{e1, e2, e3} = {Sin[rot + 2 Pi (# – 1)/3], -Cos[rot + 2 Pi (# – 1)/3]} & /@ {1, 2, 3};
{c1, c2, c3} = 1 – 3 (Sin[theta] Cos[phi – 2 Pi (# – 1)/3 – rot])^2 & /@ {1, 2, 3};
fq = Total[#[[1]] Exp[I q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3, e3}}];
{theta, phi} = {Pi/2, Pi/2};
q = {qx, qy};
fqReal = Total[#[[1]] Cos[q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3, e3}}];
fqImag = Total[#[[1]] Sin[q.#[[2]]] & /@ {{c1, e1}, {c2, e2}, {c3, e3}}];
Block[{rot = 0},
ics = Solve[
fqReal == 0 && fqImag == 0 && -Pi <= qx <= Pi && -Pi <= qy <= Pi, {qx, qy}, Reals, Method -> Reduce]
];

eqns = Equal @@@
First@Solve[{D[{fqReal, fqImag} /. {f : qx | qy :> f[rot]},
rot] == {0, 0}}, {qx'[rot], qy'[rot]}];
sols = First @ NDSolve[
{eqns, Through[{qx, qy}[0]] == ({qx, qy} /. #)},
{qx, qy},
{rot, -Pi, Pi}] & /@ ics;

NDSolve::ndsz: At rot == -0.178075, step size is effectively zero; singularity or stiff system suspected. >>
NDSolve::ndsz: At rot == 0.1780752220309632`, step size is effectively zero; singularity or stiff system suspected. >>

Some solutions go off to infinity at the singularity near Â±0.178, so I scaled down the domain of rot a bit to get a well-proportioned plot.

ParametricPlot[{qx[rot], qy[rot]} /. sols // Evaluate,
Evaluate@Flatten[{rot, 0.995 (qy[“Domain”] /. sols[[1]])}],
PlotRange -> All] /.

It appears that the set of possible positions {qx, qy} might be described as a single curve and its mirror image. The computed solutions break up the curves where the derivative of {qx, qy} with respect to rot reverses direction.

Amazing. This is very helpful, and in-fact I do need to plot a parametric plot exactly as you have. I am still astounded that there are people with such expertise on this site, willing to solve (in such helpful detail) stranger’s problems. Thank you! Maybe I can actually deliver a useful result to my supervisor for once =S
– Tom
Jan 4 ’15 at 22:26

ps there are values of ‘rot’ for which the equations really do mathematically have no solution, which is probably the reason for the singularities you found. There is a simple condition to eliminate these cases which I forgot to include in this post.
– Tom
Jan 4 ’15 at 22:28

@Michael E2: very interesting idea, that with the differential equation. I have pursued a different path considering the zeroes of fq2a = Abs[fq]^2=fqReal^2+fqImag^2. Introducing polar coordinates qx = t Cos[a], qy = t Sin[a] I then made several Plot3D’s for given t, to see the pattern zeroes, and employed Minimize f2qa with respect to rot and a for given t. I found that t>~1.35 is necessary for fq2a=0, in approximate agreement with your graph. I have not seen your critical value of rot=0.174. As Minimize Returns just a rather arbitrary Point (rot,a), the ListPlot is not a nice curve.
– Dr. Wolfgang Hintze
Jan 4 ’15 at 22:44

@Michael E2 continued: I shall edit my solution to show this in more detail. But I think you have already come very close to the final results. Also, from your graph we deduce that introducing polar coordinates for qx, qy might be not the best idea.
– Dr. Wolfgang Hintze
Jan 4 ’15 at 22:45

@Michael E2: It can be shown easily using MMA that fqReal(rot+-Pi) = fqReal(rot) and fqImag(rot+-Pi) = – fqImag(rot). Hence your chioce of the region -Pi 0)}, {qx, -d, d}, {qy, -d, d},
AxesLabel -> {qx, qy}, Axes -> True, PlotLabel -> rr, PlotPoints -> 50]
(* 150108_plot_roots_1.jpg*)

Now for rot = [Pi]/30 in the complete region g:

eps = 0.05;
d = \[Pi];
rr = \[Pi]/30;
ContourPlot[{eps == (fq2 /. rot -> rr)}, {qx, -d, d}, {qy, -d, d},
AxesLabel -> {qx, qy}, Axes -> True,
PlotLabel -> “Root positions for rot = Pi/30”, PlotPoints -> 50]

We see that the two roots z1 and z2 start to move towards each other. Also we can see that, by symmetry, is it sufficient to study the first quadrant (g1).
The next diagram, for rot = \[Pi]/18 ~= 0.174533, show the merging of z1 and z2 at a point we call z3.

eps = 0.05;
d = \[Pi];
rr = \[Pi]/18;
ContourPlot[{eps == (fq2 /. rot -> rr)}, {qx, -d, d}, {qy, -d, d},
AxesLabel -> {qx, qy}, Axes -> True,
PlotLabel -> “Root positions for rot = Pi/18”, PlotPoints -> 50]
(* 150108_plot _roots _ 3.jpg*)

Increasing rot beyond ~\[Pi]/18, there is no root left. This is the effect which Michael E2 observed with a completely different method from ours below.

Now we calculate numerically the positions z1 and z2 as a function of rot.
Our method uses Minimize[] in order to “catch” the roots.

We start “naively” with a restriction to g1:

q0 = Table[{rr,
Minimize[{fq2 /. rot -> rr, 0 <= qx < \[Pi], -0.1 <= qy < \[Pi]}, {qx, qy}]}, {rr, 0, \[Pi]/6, 0.001}]; q0p = Select[q0, #[[2, 1]] < 10^-6 &]; q0xy = {qx, qy} /. (#[[2, 2]] & /@ q0p); ListPlot-ing already give the general picture, the "root curve" qy = yq[qx], although we don't see which root is in which position of the curve ListPlot[{q0xy}, PlotRange -> {{0, 2.5}, {0, 2.5}},
PlotLabel -> “Root curve, joined picture”, AxesLabel -> {“qx”, “qy”}]
(* 150108_plot_roots_curve_1.jpg *)

The selection of a specific root (z1 or z2) is now done by restricting the region for Minimize.

For z1
(in two parts for improving accuracy)

q00 = Table[{rr,
Minimize[{fq2 /. rot -> rr, 1.3 <= qx < 1.5, -0.1 <= qy < 1.}, {qx, qy}]}, {rr, 0, 0.1, 0.001}]; q00p = Select[q00, #[[2, 1]] < 10^-6 &]; q00xy = {qx, qy} /. (#[[2, 2]] & /@ q00p); q1 = Table[{rr, Minimize[{fq2 /. rot -> rr, 0 < qx < 1.6, 0 <= qy < \[Pi]}, {qx, qy}]}, {rr, 0, \[Pi]/16, 0.001}]; q1p = Select[q1, #[[2, 1]] < 10^-6 &]; q1xy = {qx, qy} /. (#[[2, 2]] & /@ q1p); For z2 q2 = Table[{rr, Minimize[{fq2 /. rot -> rr, 1.6 < qx < \[Pi], 0 < qy < \[Pi]}, {qx, qy}]}, {rr, 0, \[Pi]/16, 0.001}]; q2p = Select[q2, #[[2, 1]] < 10^-6 &]; q2xy = {qx, qy} /. (#[[2, 2]] & /@ q2p); Before showing the complete root curve in one plot, we provide a fit to the curve Mathematica has difficulties with Fit or FindFit in this case. Therefore I did it "by eye and ruler". With the ruler I found that qx = w ~= 2.5 is a good center position for a circle for small qy. Hence the equation for qy should look like qy = Sqrt[u^2 - qx^2 - 2 w (u - qx)] /. w -> 2.5

Calculating the correcponding table qc and comparing it with the numerical curve, I found agreement after applying a factor 1.84 to the square root.

qc = Table[{qx, 1.84 Sqrt[u^2 – qx^2 – 2 w (u – qx)] /. w -> 2.5}, {qx,
u + 0.01, v, 0.01}] // Chop;

Here ist the complete plot

ListPlot[{qc, q00xy, q0xy, q1xy, q2xy}, PlotRange -> {{0, 2.5}, {0, 2.5}},
(* detailes suppressed *) AxesLabel -> {“qx”, “qy”}]
(* 150108_plot_roots_curve_02.jpg *)

This curve is to be completed by reflexions about qx- and qy-axes, independently.
The curve is in good agreement with the curve of Michael E2,
but I believe that the upper part of his curve (beyond the point of inflexion) is spurious and probably due to accuracy deficiencies in NDSolve.

In the next EDIT (#2) I shall present

(i) parametric representation of the components

{qx[rot],qy[rot]}

(ii) approximate analytic expression for these

Best regards,
Wolfgang

EDIT #2, 10.08.15

Approximate analytical treatment.

Close to rot = 0 we can derive the approximate dependence of {qx,qy} on rot analytically.

The idea is to develop the equation (letting rot -> r, for simplicity)

0 == fq2(r,qx(r),qy(r))

in a powerseries aout r = 0. As the equation must hold identically, all coefficients of the development must vanish. These conditions are then used to determine the parameters of {qx[r], qy[r]}

Let’s do it for the two roots z1 and z2 separately

For z1 (lower root)

we set, with two parameters a and b:

qx = u + a r^2
qy = b r

Skipping the details of the calculation we obtain finally

qx1[r_] = (4 ArcTan[Sqrt[3/7]])/Sqrt[3] +
2/525 r^2 (342 Sqrt[7] – 175 Sqrt[3] ArcTan[Sqrt[3/7]]);
% // N

(* Out[45]= 1.33862 + 2.77773 r^2 *)

qy1[r_] = 2/15 r (9 Sqrt[7] + 10 Sqrt[3] ArcTan[Sqrt[3/7]]);
% // N

(* Out[44]= 4.51352 r *)

Eliminating r we find the equation between qx and qy

qx1[qy_] = (4 ArcTan[Sqrt[3/7]])/Sqrt[3] + (
3 qy^2 (342 Sqrt[7] – 175 Sqrt[3] ArcTan[Sqrt[3/7]]))/(
14 (9 Sqrt[7] + 10 Sqrt[3] ArcTan[Sqrt[3/7]])^2);
% // N

(* Out[48]= 1.33862 + 0.136351 qy^2 *)

For z2 (upper root)

we set, with two parameters a1 and b1:

qx = v – a1 r
qy = h – b1 r

here

h = 2 Pi/3

Skipping the details of the calculation we obtain finally

qx2[r_] = -((2 \[Pi] r)/3) + (4 ArcTan[Sqrt[7/3]])/Sqrt[3] ;
%// N

(* Out[50]= 2.28898 – 2.0944 r *)

qy2[r_] = (2 \[Pi])/3 – 1/30 r (36 Sqrt[7] – 40 Sqrt[3] ArcTan[Sqrt[7/3]]);
%// N

(* Out[51]= 2.0944 – 0.885923 r *)

Eliminitating r we find

qy2[qx_] = (2 \[Pi])/
3 – ((9 Sqrt[7] – 10 Sqrt[3] ArcTan[Sqrt[7/3]]) (-3 qx +
4 Sqrt[3] ArcTan[Sqrt[7/3]]))/(15 \[Pi]);
% // N

(* Out[54]= 2.0944 – 0.140999 (6.86693 – 3. qx) *)

For z3 (where z1 and z2 merge for some r ~= Pi/18)

Unfortunately, I have not yet found the approach for this case, which seems to be the most interesting.

Regards,
Wolfgang

This solves my problem, thank you. Although I still have one remaining problem. If I set ‘rot’ to anything other than 0 I get the error ‘Reduce::nsmet: This system cannot be solved with the methods available to Reduce.’
– Tom
Jan 4 ’15 at 16:39

@Michael E2: I have now completed EDIT #1 (numerical treatment with Minimize) and EDIT #2 (approximate analytical treatment close to rot = 0). Please check when you find the time. I have not yet found the approximate solution close to z3 (where z1 and z2 merge) for rot = r3 (~= Pi/18 ?).
– Dr. Wolfgang Hintze
Jan 10 ’15 at 21:38

Wow, big edit! Thank you, there looks like a lot of interesting content that I will digest. It seems you’ve found the inversion symmetry and three-fold rotational symmetry of the system which is encouraging. Purely for interest, perhaps I should share that this Mathematica problem relates to finding the ‘Dirac points’ of ‘graphene-like metamaterial’ under arbitrary rotations. This work follows on from this paper currently under review: arxiv-web3.library.cornell.edu/abs/1411.7796
– Tom
Jan 12 ’15 at 11:17