7
$\begingroup$

I have the following Mathematica code:

par = {2.16023, 6.26501, 2.54649, 21.1988}; {k1, k2, k3, k4} = 
 Rationalize[par, 10^-30];

f = k1  x^2 + k2  y^2 - k3  x^3 - k4  x  y^2;

sol1 = Solve[{Grad[f, {x, y}] == 0, Det[D[f, {{x, y}, 2}]] > 0}, {x, 
    y}, Reals] // N

sol2 = Solve[{Grad[f, {x, y}] == 0, Det[D[f, {{x, y}, 2}]] < 0}, {x, 
    y}, Reals] // N

Show[Plot3D[f, {x, -0.1, 0.7}, {y, -0.25, 0.25}, Mesh -> {30}, 
  MeshFunctions -> {#3 &}, PlotTheme -> "Marketing", 
  PlotStyle -> {Cyan, Opacity[0.7]}, PlotRange -> Automatic, 
  Axes -> True, AxesLabel -> {"x", "y", "V"}, ImageSize -> 600, 
  ImageMargins -> 10, LabelStyle -> Directive[Black, 20]], 
 Graphics3D[{Red, PointSize[0.03], Point[{x, y, f} /. sol1[[1]]]}], 
 Graphics3D[{Green, PointSize[0.03], Point[{x, y, f} /. sol1[[2]]]}], 
 Graphics3D[{Blue, PointSize[0.03], Point[{x, y, f} /. sol2]}]]

which gives me the following output:

stat-pts

where red is the local minimum, green is the local maximum, and blue is the saddle points. This is one trivial example, but sometimes the expression can be complex. Is there a way in Mathematica by which I can find the stationary points directly from the plot or any other simple procedure to do that?

Improve this question
$\endgroup$
9
  • $\begingroup$ If there were a simple procedure, people would not have developed calculus. While not appreciated by many, there are good reasons why math exists $\endgroup$
    – yarchik
    Commented Jun 13 at 20:48
  • $\begingroup$ @yarchik I know that, and it's true. Actually, I have long ago seen such an answer where they extracted the stationary point from an image directly. I would also like to add calculus is definitely the basis but we also build on the great work of others to expand knowledge. $\endgroup$
    – codebpr
    Commented Jun 14 at 2:46
  • $\begingroup$ Its like you ask to drive a car on 1 wheel: It is natural to use all information about the function, like its mathematical form $f(x,y)$, to find stationary points. This I would call driving a car on 4 wheels. Still a legitimate question is how to get stationary points given a mesh representation of the function. In some cases you might miss a point, but generally it works well as t@tad demonstrates. It is like driving on 2 wheels. But I do not understand the need to discard even this information and start from just a raster image of $f(x,y)$, this is like asking for driving a car on 1 wheel. $\endgroup$
    – yarchik
    Commented Jun 14 at 4:40
  • 1
    $\begingroup$ @codebpr why can't you do this numerically without plotting it first? $\endgroup$
    – AccidentalTaylorExpansion
    Commented Jun 14 at 8:54
  • 1
    $\begingroup$ @AccidentalTaylorExpansion because I already have the plots. Well, I can go back to expression and do the numerics; that's one way, or I can step ahead and find the stationary points from the plot itself. The main motive was to learn a new way of finding that in MMA. $\endgroup$
    – codebpr
    Commented Jun 14 at 13:18

2 Answers 2

Reset to default
8
$\begingroup$

You can use one of the findAllRoots2D functions. Here's one:

(*"https://mathematica.stackexchange.com/a/97898"*)
ClearAll[findAllRoots2D];
Options[findAllRoots2D] = Join[Options[FindRoot], Options[Plot3D]];

findAllRoots2D[{f1_, f2_}, {x_, a_, b_}, {y_, c_, d_}, opts___] := 
  Module[{f1plot, f2plot}, 
   f1plot = 
    Plot3D[f1, {x, a, b}, {y, c, d}, 
     MeshFunctions -> {Function @@ {{x, y}, f1}}, Mesh -> {{0}}, 
     PlotStyle -> None, PlotRange -> All, BoundaryStyle -> None, 
     Method -> Automatic, 
     Evaluate@FilterRules[{opts}, Options[Plot3D]]];
   f2plot = 
    ListLinePlot[
     Cases[Normal@f1plot, Line[pts_] :> pts[[All, {1, 2}]], Infinity],
      MeshFunctions -> {Function @@ {{x, y}, f2}}, Mesh -> {{0}}];
   Quiet[
      Check[FindRoot[{f1 == 0, f2 == 0}, {x, #[[1]], a, 
         b}, {y, #[[2]], c, d}, 
        Evaluate@FilterRules[{opts}, Options[FindRoot]]], 
       Unevaluated@Sequence[], FindRoot::reged], FindRoot::reged] & /@
     Cases[Normal@f2plot, Point[p_] :> p, Infinity]];

cps = Union[
 findAllRoots2D[{D[f, x], D[f, y]}, {x, -0.1, 0.7}, {y, -0.25, 0.25}],
 SameTest -> 
  Function[{sol1, sol2}, 
   With[{p1 = {x, y} /. sol1, p2 = {x, y} /. sol2}, 
    Norm[p1 - p2]/(10^-16 + Max[Norm[p1], Norm[p2]]) < 10^-14]]
 ]
(*
{{x -> -1.32074*10^-25, y -> 0.},
 {x -> 0.295536, y -> -0.169578},
 {x -> 0.295536, y -> 0.169578},
 {x -> 0.565544, y -> 0.}}
*)

You can add the critical points to the plot with this:

h = D[f, {{x, y}, 2}];
colorscheme = <|{-1, -1} -> Green, {-1, 1} -> Blue, {1, 1} -> Red|>;
Graphics3D[{colorscheme[Sign@Eigenvalues[h /. #]], PointSize[0.03], 
    Point[{x, y, f} /. #]} & /@ cps]
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you so much for such an amazing answer. I would like to point out a small edit: the variable cps is not defined in the code starting with Union. But I got the result perfectly. Thanks again! $\endgroup$
    – codebpr
    Commented Jun 14 at 7:14
  • 1
    $\begingroup$ @codebpr Thanks. I added the second code later and forgot to update the first one to include setting cps. $\endgroup$
    – Michael E2
    Commented Jun 14 at 14:26
6
$\begingroup$

You can extract the points used in the 3D plot with

values = Cases[plot, GraphicsComplex[coords_, __] :> coords, \[Infinity]]

Instead of getting the points from the plot, you could create them directly, using steps small enough to resolve the curvature of your function, e.g.,

values = Flatten[Table[{{x, y}, f}, {x, -0.1, 0.7, 0.02}, {y, -0.25, 0.25, 0.02}], 1];

Then you could use those points to create an interpolating function and solve for where its gradient is zero. This is a numerical procedure so better to use FindRoot instead of Solve, assuming you have a rough idea of where the zeros are to use as initial values. This should be faster than solving symbolically with Solve when the expressions are complicated.

interp = Interpolation[values][x,y];
grad = Grad[interp, {x, y}];
det = Det[D[interp, {{x, y}, 2}]];

points = Table[
   soln = FindRoot[Thread[grad == 0], {{x, init[[1]]}, {y, init[[2]]}}];
   Append[soln, "det" -> det /. soln],
  {init, {{0.1, 0}, {0.5, 0}, {0.3, -0.2}, {0.3, 0.2}}}
]

which gives the coordinates of the points and the determinant so you can check the sign:

{{x -> 2.35223*10^-18, y -> 0., "det" -> 54.1355},
 {x -> 0.565544, y -> -2.02738*10^-18, "det" -> 49.4594}, 
 {x -> 0.295536, y -> -0.169578, "det" -> -51.6919}, 
 {x -> 0.295536, y -> 0.169578, "det" -> -51.6919}}
Improve this answer
$\endgroup$
4
  • $\begingroup$ How do you pick init, that is, the initial points ` {{0.1, 0}, {0.5, 0}, {0.3, -0.2}, {0.3, 0.2}}`, if you haven't solved the problem previously? Plot the graph, look at it, and guess? Or is there a way to do it programmatically? $\endgroup$
    – Michael E2
    Commented Jun 13 at 23:13
  • $\begingroup$ It depends on the complexity of your function. In general, finding roots of complicated nonlinear functions is difficult and you can't expect automatic procedure to always work. If you know enough about the function, e.g, based on the application, to suggest roughly where the roots are, you can use those values for 'init'. My suggestion mainly addresses speeding up the evaluation if symbolic is too slow. Not so much for how to find good initial values. Same question about the choice of range in your plot or in my Table for 'values': how do you know to pick that? $\endgroup$
    – tad
    Commented Jun 14 at 0:06
  • $\begingroup$ @tad thank you for such an informative answer. I have seen somewhere on the stack where they use MinDetect or MaxDetect for a binary image. Can such a procedure be applied to 3d plots as well? $\endgroup$
    – codebpr
    Commented Jun 14 at 2:53
  • 1
    $\begingroup$ Instead of 3D plot, you could use Density plot with grayscale gradient color function. That will give a 2D plot you can then rasterize and apply MinDetect and MaxDetect. You'd have to convert from image pixel coordinates back to the x,y values. Alternatively, Min/MaxDetect apply to an array of function values, which avoids the plot and rasterize steps. $\endgroup$
    – tad
    Commented Jun 14 at 17:45

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