12
$\begingroup$

To find all (global and local) extrema of a function in $\mathbb R^3$, I have written the following.

Example function:

n = 2.;

terrain[x_, y_] :=  2 (2 - x)^2 Exp[-(x^2) - (y + 1)^2] - 
  15 (x/5 - x^3 - y^3) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2];

fun = terrain[x, y];

plot = Plot3D[fun, {x, -n, n}, {y, -n, n}, PlotRange -> All,
              ColorFunction -> "DarkTerrain", Mesh -> False,
              PlotStyle -> Opacity@0.7]

plot of terrain function

One can observe 3 maxima and 3 minima.

NMaximize[fun, {x, y}]
{6.4547, {x -> -0.3593, y -> -0.5519}}

And

FindMaximum[fun, {x, y}]
{6.1972, {x -> -0.0529, y -> 1.2130}}

returns two of the maxima, but misses the third. My idea then was to map NMaximizeover "sufficient sectors" of the function:

 p = Flatten /@ Tuples[Partition[Range[-n, n], 2, 1], 2]
{{-2., -1., -2., -1.}, {-2., -1., -1., 0.}, ... , {1., 2., 1., 2.}}

(This algorithm was kindly provided by Kuba)

The next steps are:

max1 = NMaximize[{fun, p[[#, 1]] <= x <= p[[#, 2]], p[[#, 3]] <= y <= p[[#, 4]]},
                 {x, y}] & /@ Range@Length@p;
max2 = Chop@Partition[Cases[max1, _Real, Infinity], 3];

The result contains wrong points at the edges of the sectors, which can be deleted with

filter = # || (# /. b -> c) &[Or @@ MapThread[Equal,
         {Table[b, {n*2 + 1}], Range[-n, n]}]]
b == -2. || b == -1. || b == 0. || b == 1. || b == 2. || c == -2. ||  
c == -1. || c == 0. || c == 1. || c == 2.
max3 = DeleteCases[max2, {_, b_, c_} /; Evaluate@filter]
{{6.45471, -0.359311, -0.551929}, {6.19724, -0.0529807, 1.21301},
 {5.4426, 1.26211, -0.0152309}}

which now gives us the three maxima.

maxpoints = Graphics3D[{PointSize@0.05, Point /@ RotateLeft /@ max3}]

Repeating max1 through max3 with NMinimize finally gives this image:

density plot with extrema

Summing - up:

extrema[foo_, maxmin_, color_] :=
 Module[{res},
  res = maxmin[{foo, p[[#, 1]] <= x <= p[[#, 2]], 
       p[[#, 3]] <= y <= p[[#, 4]]}, {x, y}] & /@ Range@Length@p;
  res = Chop@Partition[Cases[res, _Real, Infinity], 3];
  res = DeleteCases[res, {a_, b_, c_} /; Evaluate@filter];
  Graphics3D[{color, PointSize@0.05, Point /@ RotateLeft /@ res}]]

Show[plot, extrema[fun, NMaximize, Black], 
 extrema[fun, NMinimize, Red], ViewPoint -> {0, 0, Infinity}]

Although my approach works, it is pretty slow (more than 2 seconds to find the extrema); and, having found it only by trial and error, I am not sure if this solution is general enough.

I would welcome any comments on how to improve this.

$\endgroup$
11
$\begingroup$
Clear["Global`*"]
n = 2.;
terrain[x_, y_] := 2 (2 - x)^2 Exp[-(x^2) - (y + 1)^2] - 
    15 (x/5 - x^3 - y^3) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2];
sol[x0_, y0_] := {x, y} /. FindRoot[
    Evaluate@{D[terrain[x, y], x] == 0, D[terrain[x, y], y] == 0}, {x,x0}, {y, y0}];
d = 0.5;
data = Table[sol[x0, y0], {x0, -n, n, d}, {y0, -n, n, d}] // Flatten[#, 1] & //
    Select[#, Function[num, Max@Abs@num < n]] & //
    DeleteDuplicates@Round[#, 10.^-6] & // Quiet;
secx[x_, y_] := Evaluate[D[terrain[x, y], {x, 2}]];
secy[x_, y_] := Evaluate[D[terrain[x, y], {y, 2}]]
secxy[x_, y_] := Evaluate[D[terrain[x, y], {x, 1}, {y, 1}]]
delta[x_, y_] := secx[x, y] secy[x, y] - secxy[x, y]^2
min = Select[data, delta @@ # > 0 && secx @@ # > 0 && secy @@ # > 0 &];
max = Select[data, delta @@ # > 0 && secx @@ # < 0 && secy @@ # < 0 &];
ContourPlot[terrain[x, y], {x, -n, n}, {y, -n, n}, Contours -> 20, 
  PlotLegends -> Automatic, ImageSize -> 300, 
  Epilog -> {Blue, PointSize[0.03], Point[min], Red, Point[max]}]

contour plot with critical points

NSolve can not solve your functions, so I can only use FindRoot to find the maxima and minima.

$\endgroup$
  • $\begingroup$ thank you so much - your solution is much faster than my attempt, it is more "mathematical", and therefore easier to read. $\endgroup$ – eldo Jun 16 '14 at 21:23
  • $\begingroup$ I have a function that is set up by a combination of different functions employing associations and switch statements, some of them use delayed evaluation (:=) and may contain vectors ({{a},{b},{c}}). All of the presented approaches produce only an empty list as a result {} (after longer or shorter computation time). Within the given range there are extrema that should be found. Do associations, switch or delayed evaluation cause those problems and how can I avoid them? Any hint (e.g. a link to an explanation) would be highly appreciated. $\endgroup$ – riddleculous Sep 28 '18 at 8:19
13
$\begingroup$

Not ideal but just for fun.

fun[a_, b_] := {x, y} /. 
  FindRoot[D[terrain[x, y], {{x, y}}] == {0, 0}, {{x, a}, {y, b}}]
h[a_, b_] := D[terrain[x, y], {{x, y}, 2}] /. {x -> a, y -> b};
pts = DeleteDuplicates[fun @@@ Tuples[Range[-2, 2, 0.5], 2]];
ptsp = Pick[pts, -2 < #[[1]] < 2 && -2 < #[[2]] < 2 & /@ pts];
col[x_, y_] := 
 If[Det[h[x, y]] < 0, {Yellow, PointSize[0.02], Point[{x, y}]}, 
  If[h[x, y][[1, 1]] > 0, {Red, PointSize[0.02], 
    Point[{x, y}]}, {Green, PointSize[0.02], Point[{x, y}]}]]
col3[x_, y_] := 
 If[Det[h[x, y]] < 0, {Yellow, PointSize[0.02], 
   Point[{x, y, terrain[x, y]}]}, 
  If[h[x, y][[1, 1]] > 0, {Red, PointSize[0.02], 
    Point[{x, y, terrain[x, y]}]}, {Green, PointSize[0.02], 
    Point[{x, y, terrain[x, y]}]}]]

Visualizing:

cp = ContourPlot[terrain[x, y], {x, -2, 2}, {y, -2, 2}, 
   Contours -> 10, Epilog -> col @@@ ptsp, 
   ColorFunction -> "DarkTerrain", ImageSize -> 300];
p3 = Show[
   Plot3D[terrain[x, y], {x, -2, 2}, {y, -2, 2}, Mesh -> False, 
    ColorFunction -> "DarkTerrain", PlotStyle -> Opacity[0.7]], 
   Graphics3D[col3 @@@ ptsp], ImageSize -> 300];
Framed@Row[{cp, p3, 
   PointLegend[{Yellow, Red, Green}, {"Saddle", "Local Min", 
     "Local Max"}]}]

enter image description here

Uses:

  • FindRoot to find critical points
  • Filtering results by DeleteDuplicatesand constraining zeros to $[-2,2]\times[-2,2]$ (to avoid 'the flatlands')
  • using second partial derivative test to classify (by color) critical points
$\endgroup$
  • 1
    $\begingroup$ +1, Your solution reminded me of a painting as an art and as an inspiration for a different approach. Thank you! Oh yeah, the painting, The Scream. $\endgroup$ – user9660 Oct 26 '15 at 8:50
  • 1
    $\begingroup$ @Lou yes, now that you comment on it, the contour plot with DarkTerrain is reminiscent of the face in 'The Scream' or the tilted head of an ET like alien! Thanks for +1 and observation :) $\endgroup$ – ubpdqn Oct 26 '15 at 8:53
  • $\begingroup$ Thanks for this new approach and for adding the saddle points, +1 $\endgroup$ – eldo Oct 26 '15 at 10:05
  • $\begingroup$ @eldo it is very similar to other answer but just thought it would be useful to exploit way to calculate Hessian etc. Thank you for +1 :) $\endgroup$ – ubpdqn Oct 26 '15 at 10:07
9
$\begingroup$

Here is my modest contribution. The idea is to use the MeshFunctions option of ContourPlot[] (as previously shown here) to extract the critical points for polishing with FindRoot[]. The Hessian is then evaluated at these points, and then tested for definiteness to identify what kind of critical points they are.

terrain[x_, y_] := 2 (2 - x)^2 Exp[-(x^2) - (y + 1)^2] -
        15 (x/5 - x^3 - y^3) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2]

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

crit = Cases[Normal[ContourPlot[dx[x, y] == 0, {x, -2, 2}, {y, -2, 2}, 
                                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}}]), ∞];

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

mini = Pick[crit, mnp]; maxi = Pick[crit, mxp]; sadl = Pick[crit, sdp];

{Legended[ContourPlot[terrain[x, y], {x, -2, 2}, {y, -2, 2}, 
                      ColorFunction -> "DarkTerrain", Contours -> 10, 
                      Epilog -> {AbsolutePointSize[6], {Cyan, Point[mini]},
                                 {Yellow, Point[sadl]}, {Orange, Point[maxi]}}], 
          PointLegend[{Cyan, Yellow, Orange}, {"Minima", "Saddles", "Maxima"}]],
 Show[Plot3D[terrain[x, y], {x, -2, 2}, {y, -2, 2}, 
             BoundaryStyle -> None, Boxed -> False, 
             ColorFunction -> "DarkTerrain", Mesh -> 10, MeshFunctions -> {#3 &}],
      Graphics3D[{{Cyan, Sphere[mini, 1/20]}, {Yellow, Sphere[sadl, 1/20]},
                  {Orange, Sphere[maxi, 1/20]}} /. {x_?NumericQ, y_?NumericQ} :>
                  {x, y, terrain[x, y]}]]} // GraphicsRow

she'll be comin' 'round the mountain...

$\endgroup$
  • $\begingroup$ I was going to do something like this, but in the process, I came up with a better answer to one of your questions $\endgroup$ – Michael E2 Oct 26 '15 at 18:07
  • $\begingroup$ Just saw it; will upvote when I can do so again. $\endgroup$ – J. M. will be back soon Oct 26 '15 at 18:09
  • $\begingroup$ @J.M. I really like this...thank you for yet another lesson :-) $\endgroup$ – ubpdqn Oct 26 '15 at 22:16
  • $\begingroup$ @ubpdqn, you're welcome! This site's great for scratching my teaching itch. :) $\endgroup$ – J. M. will be back soon Oct 26 '15 at 22:18
  • 2
    $\begingroup$ Not that important but one thing I was going to do different is Total@*Sign@*Eigenvalues /@ hl and use ±2, 0 for the Pick selection patterns for the min/max/saddle. $\endgroup$ – Michael E2 Oct 26 '15 at 23:03

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.