# Finding all maxima and minima of a function

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]]

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:

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.

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]}]

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

• thank you so much - your solution is much faster than my attempt, it is more "mathematical", and therefore easier to read.
– eldo
Jun 16, 2014 at 21:23
• 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. Sep 28, 2018 at 8:19

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"}]}]

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
• +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.
– user9660
Oct 26, 2015 at 8:50
• @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 :) Oct 26, 2015 at 8:53
• Thanks for this new approach and for adding the saddle points, +1
– eldo
Oct 26, 2015 at 10:05
• @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 :) Oct 26, 2015 at 10:07

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 *)

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]},