# 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 = saddle.nc
gmt grdmath saddle.nc DUP EXTREMA 1 EQ MUL DUP MUL 0 NAN = z.nc

gmt grd2xyz z.nc -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 saddle.nc 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");

true


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 = saddle.nc
gmt grdmath saddle.nc EXTREMA 0 NAN = b.nc
gmt grd2xyz -s b.nc
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.