Finding saddle points with grdmath

I’m trying to find saddle points with grdmath’s EXTREMA, but failing. +1/-1 both give no points (from docs: +1/-1 is saddle with max/min in x).

Here’s the code:

# make grid and extract saddle point(s)
gmt grdmath -R-10/10/-10/10 -I0.1 X X MUL Y Y MUL SUB =
gmt grdmath DUP EXTREMA 1 EQ MUL DUP MUL 0 NAN =

gmt begin saddle png
  gmt grdimage
  gmt grd2xyz -s | tee max.gmt | gmt plot
gmt end show

rm gmt.history

And the figure (just to illustrate that there is a saddle point):

Exact comparisons of floating points is never a good idea. Probably the grids has 0.9999999…

But the 1 is ‘created’ by the EXTREMA function, I thought? That’s what I’m EQ’ing for; not a value of 1 in the original grid.

EXTREMA is getting its values from the grid, which is made up of floating points, so you are comparing float == integer and 0.9999... == 1 is false

But this is just me guessing that the grid does not have a 1.0 in a single node.

But Julia tells me there is

julia> sad = grdmath("-R-10/10/-10/10 -I0.1 X X MUL Y Y MUL SUB");

julia> any(sad .== 1.0)

So EXTREMA should find saddle point…?

Yes, and it does:

gmt grdmath -R-10/10/-10/10 -I0.1 X X MUL Y Y MUL SUB =
gmt grdmath EXTREMA 0 NAN =
gmt grd2xyz -s
0 0 -1
1 Like

Thanks Paul; you’re of course right. It works.

My inexperience with grdmath, rpn (I’ve only used Casio calculators) and copy/paste of the extrema example in the docs makes my brain go full lutefisk.

… And while I’m at it; the operator D2DXY (d^2(A)/dxdy 2nd derivative); does it have a name? What is it used for? I thought it may have something to do with the $f_{xy}^2$ part of the discriminant, $f_xx \times f_yy - f_{xy}^2$, but it’s not quite it.

Not sure if it is called a cross-derivative but it is what it says. Can be useful if you need that sort of derivative of a grid to compute other stuff, e.g., elasticity and strain and that sort of thing

1 Like

Thanks Paul.