Skip to main content

plotting - 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 -> 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.



Answer



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.



Comments

Post a Comment

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more