16
$\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 -> [email protected]]

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[{[email protected], 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, [email protected], 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.

Improve this question
$\endgroup$
2

3 Answers 3

Reset to default
15
$\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.

Improve this answer
$\endgroup$
2
  • $\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, 2014 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, 2018 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
Improve this answer
$\endgroup$
4
  • 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, 2015 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, 2015 at 8:53
  • $\begingroup$ Thanks for this new approach and for adding the saddle points, +1 $\endgroup$
    – eldo
    Oct 26, 2015 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, 2015 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...

Improve this answer
$\endgroup$
11
  • $\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, 2015 at 18:07
  • $\begingroup$ Just saw it; will upvote when I can do so again. $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 26, 2015 at 18:09
  • $\begingroup$ @J.M. I really like this...thank you for yet another lesson :-) $\endgroup$
    – ubpdqn
    Oct 26, 2015 at 22:16
  • $\begingroup$ @ubpdqn, you're welcome! This site's great for scratching my teaching itch. :) $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 26, 2015 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, 2015 at 23:03

Your Answer

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

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