3
$\begingroup$

In short: I have a 2-D function that gives this plot, where the orange dots are maxima: contour plot surface plot

There are 4 maxima in this plot.

Now, I plot the function again in a smaller region, close to these maxima, resulting in this this:

So, one of the maxima turned out to be fake.

How can I change the main code to better extremise the function (when it is on a larger area, like the first picture) and not return false maxima?

Thanks so much.

Code:

k = 2 π;
k1 = {k, 0};
k2 = FullSimplify[RotationMatrix[(2 π)/5]. k1];
k3 = FullSimplify[RotationMatrix[(2 π)/5]. k2];
k4 = FullSimplify[RotationMatrix[(2 π)/5] .k3];
k5 = FullSimplify[RotationMatrix[(2 π)/5] .k4];
E1[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k1.{x, y} - α1]
E2[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k2.{x, y} - α1]
E3[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k3.{x, y} - α1]
E4[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k4.{x, y} - α1]
E5[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k5.{x, y} - α1]
Eges[x_, y_] := 
 E1[x, y, 1, 0] + E2[x, y, 1, 0] + E3[x, y, 1, 0] + E4[x, y, 1, 0] + 
  E5[x, y, 1, 0]
IIges[x_, y_] := 
  Conjugate[Eges[x, y]] Eges[x, y] // ComplexExpand // 
   Simplify; (* can't differentiate Conjugate, had to change this *)
xm = 4; 

{dx[x_, y_], dy[x_, y_]} = 
  D[IIges[x, y], {{x, y}}]; (* calculating partial derivatives *)
hes[x_, y_] = 
  D[IIges[x, y], {{x, y}, 
    2}]; (* Hessian matrix: Second derivative matrix *)

crit = Cases[
   Normal[ContourPlot[dx[x, y] == 0, {x, -xm, xm}, {y, -xm, xm}, 
     ContourStyle -> None, Mesh -> {{0}}, 
     MeshFunctions -> Function[{x, y, z}, dy[x, y]]]], 
   Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /. 
      FindRoot[{dx[\[FormalX], \[FormalY]], 
        dy[\[FormalX], \[FormalY]]}, {{\[FormalX], x0}, {\[FormalY], 
         y0}}]),  ∞];

(* Identify points as minima, maxima or saddle points *)
hl = hes @@@ crit;
mnp = PositiveDefiniteMatrixQ /@ hl; (*pick minima*)
mxp = PositiveDefiniteMatrixQ /@ (-hl); (*pick maxima*)
sdp = Thread[mnp ⊽ mxp];(*saddle points are leftovers*)

(* Finding coordinates of critical points *)
mini = Pick[crit, mnp];
maxi = Pick[crit, mxp];
sadl = Pick[crit, sdp];

(* Plot stuff!
I changed the Image Size because it was too small *)
{Legended[
  ContourPlot[IIges[x, y], {x, -xm, xm}, {y, -xm, xm}, 
   ImageSize -> Scaled[0.3], ColorFunction -> "DarkTerrain", 
   Contours -> 15, 
   Epilog -> {AbsolutePointSize[6], {Cyan, Point[mini]}, {Orange, 
      Point[maxi]}}], 
  PointLegend[{Cyan, Orange}, {"Minima", "Maxima"}]], 
 Show[Plot3D[IIges[x, y], {x, -xm, xm}, {y, -xm, xm}, 
   BoundaryStyle -> None, Boxed -> False, 
   ColorFunction -> "DarkTerrain", Mesh -> 15, 
   MeshFunctions -> {#3 &}], 
  Graphics3D[{{Cyan, Sphere[mini, 1/20]}, {Orange, 
      Sphere[maxi, 1/20]}} /. {x_?NumericQ, y_?NumericQ} :> {x, y, 
      IIges[x, y]}], ImageSize -> Scaled[0.43]]}

crit2 = Cases[
   Normal[ContourPlot[dx[x, y] == 0, {x, 0.5, 2}, {y, -2.5, -1}, 
     ContourStyle -> None, Mesh -> {{0}}, 
     MeshFunctions -> Function[{x, y, z}, dy[x, y]]]], 
   Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /. 
      FindRoot[{dx[\[FormalX], \[FormalY]], 
        dy[\[FormalX], \[FormalY]]}, {{\[FormalX], x0}, {\[FormalY], 
         y0}}]), ∞];
hl2 = hes @@@ crit2;
mnp2 = PositiveDefiniteMatrixQ /@ hl2; (*pick minima*)
mxp2 = PositiveDefiniteMatrixQ /@ (-hl2); (*pick maxima*)
sdp2 = Thread[mnp2 ⊽ mxp2];(*saddle points are leftovers*)

mini2 = Pick[crit2, mnp2];
maxi2 = Pick[crit2, mxp2];
sadl2 = Pick[crit2, sdp2];

{Legended[
  ContourPlot[IIges[x, y], {x, 0.5, 2}, {y, -2.5, -1}, 
   ImageSize -> Scaled[0.3], ColorFunction -> "DarkTerrain", 
   Contours -> 15, 
   Epilog -> {AbsolutePointSize[6], {Cyan, Point[mini2]}, {Orange, 
      Point[maxi2]}}], 
  PointLegend[{Cyan, Orange}, {"Minima", "Maxima"}]], 
 Show[Plot3D[IIges[x, y], {x, 0.5, 2}, {y, -2.5, -1}, 
   BoundaryStyle -> None, Boxed -> False, 
   ColorFunction -> "DarkTerrain", Mesh -> 15, 
   MeshFunctions -> {#3 &}], 
  Graphics3D[{{Cyan, Sphere[mini2, 1/20]}, {Orange, 
      Sphere[maxi2, 1/20]}} /. {x_?NumericQ, y_?NumericQ} :> {x, y, 
      IIges[x, y]}], ImageSize -> Scaled[0.43]]}
| improve this question | |
$\endgroup$
7
$\begingroup$

I don't quite feel like writing out the full solution. Nevertheless, I see that you're trying your best to adapt one of my previous solutions, and I will at least show a way for you to work a little bit smarter.

One of the tenets of numerical computing is the one of exploiting the structure of the problem; what this says is that you look if there is any structure like symmetry that can be used to reduce computational effort, improve numerical stability, or even both.

In your particular case, I previously gave you a hint: Conjugate[Eges[x, y]] == Eges[-x, -y]. Apart from the obvious advantage of not needing to differentiate Conjugate[] anymore, it's a clue that your objective function has twofold symmetry.

Another exploitable fact is in how your function was built: one can see that the function ought to have five-fold symmetry. Put those two facts together, and we find that the objective function is ten-fold ($\mathbf{5\times2}$) symmetric.

How to exploit this, you ask? This means we can concentrate all our computational efforts within a wedge of the plane, and later reproduce all the solutions through suitable rotations.

Thankfully, ContourPlot[] supports a RegionFunction option, which allows you to restrict effort to just a slice of the whole pie.

Here is how you can start. Note that I have simplified the generation of the objective function:

Clear[x, y];
k1 = {2 π, 0}; rots = RootReduce[NestList[RotationTransform[2 π/5], k1, 4]];
Eges[x_, y_, I1_: 1, a1_: 0] = Sqrt[I1] Total[Exp[I rots.{x, y} - a1]];
IIges[x_, y_] = Simplify[ComplexExpand[Re[Eges[-x, -y] Eges[x, y]]]];

{dx[x_, y_], dy[x_, y_]} = D[IIges[x, y], {{x, y}}];
hes[x_, y_] = D[IIges[x, y], {{x, y}, 2}];

(* note how RegionFunction "slices the pie" *)
crit = Cases[Normal[ContourPlot[dx[x, y] == 0, {x, -4 Sqrt[2], 4 Sqrt[2]},
                                {y, -4 Sqrt[2], 4 Sqrt[2]}, 
                                ContourStyle -> None, Mesh -> {{0}},
                                MeshFunctions -> Function[{x, y, z}, dy[x, y]],
                                PlotPoints -> 75, RegionFunction ->
                                Function[{x, y, z}, 0 < Arg[x + I y] < π/5 &&
                                         Sqrt[x^2 + y^2] < 4 Sqrt[2]]]], 
             Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /. 
             FindRoot[{dx[\[FormalX], \[FormalY]], dy[\[FormalX], \[FormalY]]},
                      {{\[FormalX], x0}, {\[FormalY], y0}}]), ∞];

hl = hes @@@ crit;
mnp = PositiveDefiniteMatrixQ /@ hl; (* pick minima *)
mxp = PositiveDefiniteMatrixQ /@ (-hl); (* pick maxima *)
mini = Pick[crit, mnp]; maxi = Pick[crit, mxp];

Visualize the solutions in the chosen wedge:

ContourPlot[IIges[x, y], {x, 0, 4 Sqrt[2]}, {y, 0, 4 Sqrt[2] Sin[π/5]}, 
            AspectRatio -> Automatic, ColorFunction -> "ThermometerColors", 
            Epilog -> {AbsolutePointSize[5],
                       {Orange, Point[mini]}, {Green, Point[maxi]}},
            PlotPoints -> 45,
            RegionFunction -> Function[{x, y, z}, 0 < Arg[x + I y] < π/5 &&
                                       Sqrt[x^2 + y^2] < 4 Sqrt[2]]]

a slice of the pie, with some toppings

A perceptive eye will notice two things: the origin was missed even though it was a critical point (not too hard to fix), and some of the points in one edge of the "slice" were also captured in the other edge. For that second caveat, this means that when you do the rotation later to reproduce all the critical points, you will have to strip out some duplicates. Still, not having to look at too many critical points is a big thing, and we are lucky that symmetry allows this simplification.

| improve this answer | |
$\endgroup$
  • $\begingroup$ (+1) You have generated a segment of a very nice snowflake-like texture! Festive mood! :) I would appreciate if you post the complete figure. $\endgroup$ – Alexey Popkov Nov 13 '15 at 10:23
  • $\begingroup$ In due time, @Alexey. In due time. :) $\endgroup$ – J. M.'s technical difficulties♦ Nov 13 '15 at 10:33
  • $\begingroup$ Thanks, but I did not want to assume any symmetry a priori, I would have rather wanted to see it from the plot. Thanks though! $\endgroup$ – SuperCiocia Nov 13 '15 at 12:31
  • $\begingroup$ @Super, "did not want to assume any symmetry" - well, the symmetry only was determined by looking at the properties of the function; no assumptions apart from $x,y$ being real was used. The plot merely confirms the ten-fold symmetry. In any event, it would be wasteful not to try doing some analytical legwork before plunging into numerics. $\endgroup$ – J. M.'s technical difficulties♦ Nov 13 '15 at 12:54
  • $\begingroup$ My interest actually lied in a larger plot, xm = 8. If you plot that, little pentagons arise, but sometimes they are "missing a point". In the small plot (xm=4), these 'broken pentagons' look like the figure circled in red in the first picture. Since I was not sure whether or not to expect these pentagons, I focussed on a small region (this question) and tried to see whether the 'extra points' that I saw actually belonged to rhombs or pentagons. $\endgroup$ – SuperCiocia Nov 13 '15 at 14:08
8
$\begingroup$

The extraneous points result from inadequate resolution in the ContourPlot used to provide initial estimates for crit. A simple but effective approach is to increase PlotPoints.

ContourPlot[dx[x, y] == 0, {x, -xm, xm}, {y, -xm, xm}, PlotPoints -> 50,
    ContourStyle -> None, Mesh -> {{0}}, MeshFunctions -> Function[{x, y, z}, dy[x, y]]]

critical points

This change reduces Length[crit] from 895 to 853 and also eliminates numerous FindRoot::lstol: error messages.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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