0
$\begingroup$

I am trying to find all the local minima of a function, in a specified region of space. Unfortunately this is taking way too long for large system sizes, and I am wondering whether it's just because my approach is unnecessarily computationally heavy.

Preliminary variables and names:

\[Mu]m = 10^-6;
nm = 10^-9; (* m *)
xm = {50, 100, 500, 1500}; (* window size in units of \[Lambda] *) 
\[Lambda]Lattice = 725 nm;
kLattice = (2 \[Pi])/\[Lambda]Lattice;    
kdim = kLattice*1/\[Lambda]Lattice;
depthK = 10; 

This is the function I need to find the minima of:

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

IdealPotential[x_, y_] = 
 IdealPotential[depthK, kdim, x, y];

IdealPotentialData = 
 Table[FindExtrema[IdealPotential, xm[[n]], 250], {n, 1, Length[xm]}];

{IdealMini, IdealMinima, IdealcurvMini, IdealMaxi, IdealMaxima, 
  IdealcurvMaxi} = Transpose[IdealPotentialData ];

calculated via the method FindExtrema, defined as such:

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]]]], 
    Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /. 
       FindRoot[{dx[\[FormalX], \[FormalY]], 
         dy[\[FormalX], \[FormalY]]}, {{\[FormalX], x0}, {\[FormalY], 
          y0}}]), \[Infinity]];
  hl = hes @@@ crit;
  mnp = PositiveDefiniteMatrixQ /@ hl; 
  mxp = PositiveDefiniteMatrixQ /@ (-hl); 
  sdp = Thread[mnp \[Nor] mxp];
  mini = Pick[crit, mnp];
  maxi = Pick[crit, mxp];
  sadl = Pick[crit, sdp];
  {mini, potential @@@ mini, hes @@@ mini, maxi, potential @@@ maxi, 
   hes @@@ maxi}
  ]

-

The PlotPoints function in particular is essnetial for finding all the minima but adds tremendously to computational time.

Ideas

  • My first idea would be to only look for minima in a certain region, since my potential is symmetric around the origin. I am currently looking into that...

Does anyone have any other suggestion on how to speed this up, or on a different method? Thanks. Maybe there's a way to cast this problem so as to use NSolve?

$\endgroup$
  • 1
    $\begingroup$ Look at the result of Plot3D[IdealPotential[x, y], {x, -2, 2}, {y, -2, 2}]. $\endgroup$ – user64494 Nov 27 '17 at 19:32
  • $\begingroup$ Yes, which is why I increase the PlotPoints, but then it takes ages. $\endgroup$ – SuperCiocia Nov 27 '17 at 21:52
  • 1
    $\begingroup$ Did you try FindMinimum and NMinimize? Which options did you try? $\endgroup$ – M. Stern Nov 27 '17 at 23:32
  • $\begingroup$ From what I could tell, the distance between the local peaks of IdealPotential[x,y] is about 10^-13. That means that for the region {x, y} ∈ [0, 1]×[0, 1] you're looking for about 10^26 points... $\endgroup$ – aardvark2012 Nov 28 '17 at 10:33

Your Answer

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

Browse other questions tagged or ask your own question.