Finding saddle points with grdmath

Q&A
1

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 begin saddle png
  gmt grdimage saddle.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):

saddle
saddle1772×1772 87.1 KB

2

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

3

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.

4

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.

5

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)
true
6

So EXTREMA should find saddle point…?

7

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
8

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.

9

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
10

Thanks Paul.