0
$\begingroup$

I am having a problem similar to the previous question here. Solve either takes forever (I assume infinity) to run or returns an error. I have tried various combinations of using Round, Chop, N and Real to neaten up output, but I can't quite pinpoint where the problem is. In my previous question it was a mixture of exact and inexact coefficients causing Solve to run indefinitely, but here I am not sure. Normally I would post an analogous problem to make it easier for those reading but as I'm not sure even roughly where the issue is this time, I'm afraid I will post the full problem.

(*Nearest neighbour vectors*)
d1 = {0, -1} - {rx, ry};
d2 = {Sqrt[3]/2, 1/2} - {rx, ry};
d3 = {-Sqrt[3]/2, 1/2} - {rx, ry};

(*The coefficients*)
a = Norm[d1]^-3; b = Norm[d2]^-3; c = Norm[d3]^-3;

(*modulus of fq*)
modfq[k_] := 
  Sqrt[a^2 + b^2 + c^2 + 2 a b Cos[k.(d1 - d2)] + 
    2 b c Cos[k.(d2 - d3)] + 2 a c Cos[k.(d1 - d3)]];

(*Differential of fq w.r.t x and y to find turning points and then find minimum*)
findMin[p_] := (rx = p[[1]]; ry = p[[2]];
  Min[modfq /@ (Re[{x, y} /. 
       Solve[a b (d2 - d1)[[1]] Sin[{x, y}.(d1 - d2)] + 
          b c (d3 - d2)[[1]] Sin[{x, y}.(d2 - d3)] + 
          a c (d3 - d1)[[1]] Sin[{x, y}.(d1 - d3)] == 
          a b (d2 - d1)[[2]] Sin[{x, y}.(d1 - d2)] + 
          b c (d3 - d2)[[2]] Sin[{x, y}.(d2 - d3)] + 
          a c (d3 - d1)[[2]] Sin[{x, y}.(d1 - d3)] == 0, {x, y}]])])

(*Array of positions*)
res = 2; (*<--Ultimately this will be a value of a few thousand*)
xArray = Range[-Sqrt[3]/2, Sqrt[3]/2, Sqrt[3]/res];
yArray = Range[3/4, -3/4, -3/(2 res)];
g[p_] := Transpose@{xArray, ConstantArray[p, Length[xArray]]};
rArray = g /@ yArray ;

(*Map findMin onto rArray*)
diagram = Map[findMin, rArray, {2}];
ArrayPlot[diagram, ColorFunction -> "Rainbow"]

In this form it runs indefinitely, and attempts at rounding, choping etc. just don't seem to be working for me. Any suggestions on what the problem is?

Many thanks for any help, and apologies for the big code dump.

$\endgroup$
10
  • $\begingroup$ Solve::nsmet : "This system cannot be solved with the methods available to Solve." You can't expect magic to happen, if when using Mathematica. Some equations can't be analytically solved. $\endgroup$
    – Nasser
    Jul 4, 2014 at 12:31
  • $\begingroup$ Have you unit tested each of the functions you defined in your code, so that you know they each work? $\endgroup$
    – m_goldberg
    Jul 4, 2014 at 12:32
  • $\begingroup$ I've been testing as I can go and only ran into errors in the few final steps :\ $\endgroup$
    – Tom
    Jul 4, 2014 at 12:38
  • $\begingroup$ @m_goldberg that is what I did first. Broke that clooped up code into separate cells to see if each work. The problem is Solve can't solve the equations even with Real option. I do not know why people do not break code into separate pieces and run each in sequence any more like we used to do in the good old days. Wrapping functions inside other functions might make things look cool, but it makes debugging very hard. $\endgroup$
    – Nasser
    Jul 4, 2014 at 12:40
  • $\begingroup$ Huh, I did test the Solve with some input whilst writing and it ran and put out an answer I was expecting.. I don't know.. $\endgroup$
    – Tom
    Jul 4, 2014 at 12:43

1 Answer 1

2
$\begingroup$
Clear["Global`*"];
f[rx_, ry_] := Module[{d1, d2, d3, a, b, c, sol1, sol2, data},
    {d1, d2, d3} = # - {rx, ry} & /@ {{0, -1}, {Sqrt[3]/2, 1/2}, {-Sqrt[3]/2, 1/2}};
    {a, b, c} = Norm[#]^-3 & /@ {d1, d2, d3};
    sol1 = NSolve[a b (d2 - d1) Sin[{x, y}.(d1 - d2)] + 
       b c (d3 - d2) Sin[{x, y}.(d2 - d3)] + 
       a c (d3 - d1) Sin[{x, y}.(d1 - d3)] == {0, 0}, {x, y}][[All, All, 2]];
    sol2 = Re@Flatten[Table[sol1 /. {C[1] -> i, C[2] -> j}, {i, -1, 1}, {j, -1, 1}], 2];
    data = Sqrt[a^2 + b^2 + c^2 + 2 a b Cos[#.(d1 - d2)] + 
       2 b c Cos[#.(d2 - d3)] + 2 a c Cos[#.(d1 - d3)]] & /@ sol2;
    Min@Chop[data, 10^-7]
]
res = 20;
diagram = Monitor[Table[f[rx, ry], 
    {rx, -Sqrt[3]/2, Sqrt[3]/2, Sqrt[3]/res}, {ry, 3/4, -3/4, -3/(2 res)}], 
    N[(rx + Sqrt[3]/2)/Sqrt[3]]];
ArrayPlot[diagram, ColorFunction -> "Rainbow", ColorFunctionScaling -> False]

enter image description here

$\endgroup$
1
  • $\begingroup$ Nice. That looks like what I was expecting :) Sure seems to take a long time to run though =\ I might have to approach the physical problem a different way $\endgroup$
    – Tom
    Jul 4, 2014 at 19:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.