Let $f : D\subseteq\mathbb{R}^2 \longrightarrow \mathbb{R}^2$ be differentiable. Denote by $$\mathcal{C}=\{x\in D \,|\, \det \mathfrak{J}_f(x)=0 \}$$ the set of all critical points, where the Jacobian matrix of $f$ is singular.
Assume that $f$ is known only numerically. Is there a nice numerical method to get those points where $\det \mathfrak{J}_f \approx 0$, in particular, when the calculation of $\mathfrak{J}_f$ is elaborating?
Thanks in advance for some tips, literature, etc.!
Idea 1: Let $p \in D$. From $$f(x) - f(p) = \mathfrak{J}_f(p)(x- p)^t + o(\|x -p\|) \quad (\text{as } x \to p)$$ one can derive for $\epsilon >0$ that $$\mathfrak{J}_f(p)_{i,j} \approx\frac{f_i(p+\epsilon e_j)-f_i(p)}{\epsilon}, \ i,j\in \{1,2\}.$$ This allows to get $\det \mathfrak{J}_f$ approximately. The 'Marching squares' algorithm can then be used to find the contour where it vanishes. Unfortunately, this idea is brute-force and might require a fairly dense grid.
Idea 2: Sample points $x_1, \dots, x_n$ uniformly at random in $D$. The regions where the pattern $$\{ r(x_i) \,|\, i=1, \dots,n \}$$ is dense are regions where $|\det \mathfrak{J}_f(x_i)|$ is low due to 'Change of Variables Theorem'. Unfortunately, this might be costly as well.
Consider a point $x\in D$ and an infinitesimal rectangular triangle $S \subset D$ such that $x\in S$. The triangle is transformed by $f$ to some curved region $f(S)$. In contrast, the best linear approximation $\mathfrak{J}_f(x)$ of $f$ around $x$ maps $S$ again to a triangle $P$ with $$\text{signarea}(P) = \det \mathfrak{J}_f(x) \cdot \text{signarea}(S),$$ where $\text{signarea}$ denotes the signed area. Instead of using $P$ one might use the triangle $Q$ given by the $f$-transformed vertices of $S$. It is $$\det \mathfrak{J}_f(x) = \frac{\text{signarea}(P)}{\text{signarea}(S)}\approx \frac{\text{signarea}(Q)}{\text{signarea}(S)}.$$ If $S=\Delta(x_1,x_2,x_3)$ the signed area can be calculated using the 'Shoelace formula' by $$\text{signarea}(Q)=\frac{1}{2}\left(f_1(x_1) f_2(x_2) + f_1(x_2) f_2(x_3) + f_1(x_3) f_2(x_1) - f_1(x_2) f_2(x_1) - f_1(x_3) f_2(x_2) - f_1(x_1) f_2(x_3) \right).$$ For instance using $x_1=x+\epsilon e_1, x_2=x+\epsilon e_2, x_3=x$ leads to $$\det \mathfrak{J}_f(x) \approx F_\epsilon(x)=\frac{1}{\epsilon^2}\left(f_1(x+\epsilon e_1) f_2(x+\epsilon e_2) + f_1(x+\epsilon e_2) f_2(x) + f_1(x) f_2(x+\epsilon e_1) - f_1(x+\epsilon e_2) f_2(x+\epsilon e_1) - f_1(x) f_2(x+\epsilon e_2) - f_1(x+\epsilon e_1) f_2(x) \right).\tag{1}$$
Now a multivariate root finding algorithm can be applied to $F_\epsilon$.