# Why does my code find minima outside the square that constrains my search?

I have written a script to find the minima of a function, and then manipulate them. I have a physical expectation on what the result should look like, but I am not seeing this and I am wondering whether numerical errors are part of the problem.

I am finding the minima of a function RealPotential in the xy-plane, meaning their locations and their values. I am restricting the search to a Abs[x] <= xm/2 and Abs[y] <= xm/2, which is the window size.

I then calculate the values of another function, IdealPotential, at the same minima and plot the difference between these two sets of minima. This difference should go rouglhly as a Gaussian, with some noise around.

This difference between the two sets of minima I get is this: The plot above is obtained with the script running for a window size xm of 100, so it should only find minima between -50 and 50. But some points also appear around 200. This is weird because the function has an eightfold symmetry, which is observed by these points as well. Also, the points marked by the red arrow do not follow a Gaussian distribution.

The plot above was obtained withPlotRange -> Full (in the last line of code presented at the end of the question). When I restrict the range with PlotRange -> {{-50, 50}, {-50, 50}, {-0.1, 0.1}}, I see something that makes sense and that I expect: • How is it possible for the script to find minima outside the specified region of interest?

• Is there a way to quantify numerical errors responsible for the werid points not-Gaussian following that I see?

### Code

My method to find the minima is based on on this.

Functions to find minima of RealPotential, and calculate IdealPotential values at the same minima:

FindExtrema[potential_, windowsize_, points_] :=
Module[ {dx, dy, hl, x, y, hes, crit, mnp, mxp, sdp, mini, maxi,
sadl},
{dx[x_, y_], dy[x_, y_]} = D[potential[x, y], {{x, y}}];
hes[x_, y_] = D[potential[x, y], {{x, y}, 2}];
crit = Cases[
Normal[ContourPlot[
dx[x, y] == 0, {x, -windowsize/2, windowsize/
2}, {y, -(windowsize/2), windowsize/2}, PlotPoints -> points,
ContourStyle -> None, Mesh -> {{0}},
MeshFunctions -> Function[{x, y, z}, dy[x, y]],
RegionFunction ->
Function[{x, y, z},
0 < Arg[x + I y] < π/4 &&
Sqrt[x^2 + y^2] < windowsize/2]]],
Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /.
FindRoot[{dx[\[FormalX], \[FormalY]],
dy[\[FormalX], \[FormalY]]}, {{\[FormalX], x0}, {\[FormalY],
y0}}]), ∞];
hl = hes @@@ crit;
mnp = PositiveDefiniteMatrixQ /@ hl;
mxp = PositiveDefiniteMatrixQ /@ (-hl);
sdp = Thread[mnp ⊽ mxp];
mini = Pick[crit, mnp];
maxi = Pick[crit, mxp];
sadl = Pick[crit, sdp];
{mini, potential @@@ mini, hes @@@ mini, maxi, potential @@@ maxi,
hes @@@ maxi}
]

FindIdealValues[idealpotential_, mini_] := Module[ {x, y, hes},
hes[x_, y_] = D[idealpotential[x, y], {{x, y}, 2}];
{idealpotential @@@ mini, hes @@@ mini}
]


Data just needed to set up the scene:

μm = 10^-6;
nm = 10^-9; (* m *)
MHz = 10^6;
umass = 1.67*10^-27; (* kg *)
mRb = 87* umass;
mK = 39.5*umass;
h = 6.626070040*10^-34;
ℏ = 1.0546*10^-34;
c = 299792458;

λRbD1 = 794.98 nm;
λRbD2 = 780.24 nm;
ΓRbD1 = 2 π *5.746 MHz;
ΓRbD2 = 2 π *6.065 MHz;
ωRbD1 = (2 π c)/λRbD1;
ωRbD2 = (2 π c)/λRbD2;

λKD1 = 770.11 nm;
λKD2 = 766.7 nm;
ΓKD1 = 2 π*5.956 MHz;
ΓKD2 = 2 π*6.035 MHz;
ωKD1 = (2 π c)/λKD1;
ωKD2 = (2 π c)/λKD2;

αRb = -((π c^2)/(
2 ωRbD1^3))*ΓRbD1*(1/(ωRbD1 - \
ωLattice) + 1/(ωRbD1 + ωLattice)) - (π c^2)/(
2 ωRbD2^3)*ΓRbD2*(2/(ωRbD2 - \
ωLattice) + 2/(ωRbD2 + ωLattice));
αK = -((π c^2)/(
2 ωKD1^3))*ΓKD1*(1/(ωKD1 - \
ωLattice) + 1/(ωKD1 + ωLattice)) - (π c^2)/(
2 ωKD2^3)*ΓKD2*(2/(ωKD2 - \
ωLattice) + 2/(ωKD2 + ωLattice));

λLattice = 725 nm;
kLattice = (2 π)/λLattice;
ωLattice = (2 π c)/λLattice;
whLattice = 200 μm;
wvLattice = 80 μm;

ErRb = (ℏ^2*kLattice^2)/(2*mRb);(* Rb recoil energy*)
ErK = (ℏ^2*kLattice^2)/(2*mK); (* K recoil energy*)

V0[depth_, Er_] := depth*Er;
power[α_, Er_, depth_, wh_, wv_] :=
V0[depth, Er]/(8 α) π wh wv;

depthK = 10; (* depth of 1D lattice in recoil energies for Potassium *)

powerK = power[αK, ErK, depthK, whLattice, wvLattice] (* W *)

depthRb = (2 powerK)/(π whLattice wvLattice) *(4 αRb)/ErRb


Rescaling to dimensionless units (so as to not use microns but integers) and setting system syze:

xm = {100}; (* window size in units of λ *)
plotpoints = {100};
kdim = kLattice*λLattice;
factor = 1/Sqrt[2 π^2];
whLatticedim = whLattice/λLattice;
wvLatticedim = wvLattice/λLattice;


RealPotential definition:

RealPotential[V0_, k_, wh_, x_, y_] :=
V0*(Sin[k x]^2 *Exp[-((2 y^2)/wh^2)] +
Sin[k y]^2*Exp[-((2 x^2)/wh^2)] +
Sin[k/Sqrt (x + y)]^2*Exp[-((x - y)^2/wh^2)] +
Sin[k/Sqrt (x - y)]^2*Exp[-((x + y)^2/wh^2)]);
RealPotential[x_, y_] =
RealPotential[depthK, kdim, whLatticedim, x, y];

RealPotentialData =
Table[FindExtrema[RealPotential, xm[[n]], plotpoints[[n]]], {n, 1,
Length[xm]}];

{RealMini, RealMinima, RealcurvMini, RealMaxi, RealMaxima,
RealcurvMaxi} = Transpose[RealPotentialData ];

RealGSEnergy =
Table[1/2*
factor*(Sqrt[Eigenvalues[RealcurvMini[[n]][[m]]][]] + Sqrt[
Eigenvalues[RealcurvMini[[n]][[m]]][]]), {n, 1,
Length[xm]}, {m, 1, Length[RealcurvMini[[n]]]}]; (*+RealMinima ? *)

RealMiniRot =
Table[Join @@
NestList[RotationTransform[π/4], RealMini[[n]], 7], {n, 1,
Length[RealMini]}];


Plots (not essential):

 Table[
{Legended[
ContourPlot[
RealPotential[x, y], {x, -(xm[[s]]/0.4), xm[[s]]/
0.4}, {y, -xm[[s]]/0.4, xm[[s]]/0.4}, ImageSize -> Scaled[0.3],
ColorFunction -> "ThermometerColors", Contours -> 15,
Epilog -> {AbsolutePointSize, {Opacity, Red,
Point[RealMiniRot[[s]]]}, {Red, Point[RealMiniRot[[s]]]}}],
PointLegend[{Red}, {"Minima"}]]}, {s, 1, Length[xm]}]


Finding the values of the IdealPotential:

IdealPotential[V0_, k_, x_, y_] :=
V0*(Sin[k x]^2  + Sin[k y]^2 + Sin[k/Sqrt (x + y)]^2 +
Sin[k/Sqrt (x - y)]^2);
IdealPotential[x_, y_] = IdealPotential[depthK, kdim, x, y];

IdealPotentialData =
Table[FindIdealValues[IdealPotential, RealMini[[n]]], {n, 1,
Length[xm]}];

{IdealMinima, IdealcurvMini} = Transpose[IdealPotentialData ];

IdealMinimaGS =
Table[1/2*
factor*(Sqrt[Eigenvalues[IdealcurvMini[[n]][[m]]][]] + Sqrt[
Eigenvalues[IdealcurvMini[[n]][[m]]][]]), {n, 1,
Length[xm]}, {m, 1, Length[RealcurvMini[[n]]]}];


Plots:

DifferenceGSEnergy =
Table[(RealGSEnergy[[n]][[m]] - IdealMinimaGS[[n]][[m]]), {n, 1,
Length[xm]}, {m, 1, Length[RealcurvMini[[n]]]}];

DifferenceXYE =
Table[Flatten[{RealMini[[n]][[m]],
DifferenceGSEnergy[[n]][[m]]}], {n, 1, Length[xm]}, {m, 1,
Length[RealcurvMini[[n]]]}]

DifferenceXYERot =
Table[Join @@ (Flatten /@
Thread[{NestList[RotationTransform[π/4], {#, #2},
7], #3}] & @@@ DifferenceXYE[[n]]), {n, 1,
Length[DifferenceXYE]}];

Table[ListPointPlot3D[DifferenceXYERot[], PlotRange -> Full,
ImageSize -> 500], {j, 1, Length[DifferenceXYERot]}]

• AFAIUI (I have no time to go into details of your code.), you find critical points, solving a system of partial derivatives equated to zero by FindRoot command.This may cause your problem. Did you try the direct approach, having made use of NMInimize, FindMaximum etc? – user64494 May 15 '18 at 5:27
• The thing is all the Mathematica built in functions always aim to finding the global minimum/maximum, whereas I would like to find ALL the local extrema. – SuperCiocia May 15 '18 at 10:04
• As a first step, you may split the domain into small domains and find global extrema on these domains. – user64494 May 15 '18 at 12:30
• One more suggestion is to filter critical points outside the constraints. – user64494 May 16 '18 at 7:44