# Continuation of strictly monotone function on a plane

I have a computational problem where I'm given the values of a function $f(x,y)$ sampled on 2D regular grid. The function is given in a compact domain of a grid and it is strictly monotone w.r.t. $x$ and $y$ in the domain. It should be continued to the whole plane, with the only requirement of being continuous and strictly monotone w.r.t. every argument.

In a special case, when the domain is a square, I can construct an explicit formula for such continuation. So what I miss is a formula or an algorithm to continue the function to the bounding box of the domain, supporting the monotonic property above.

As a brute force algorithm I can imagine representing monotony conditions for all points on the grid (outside $A$, in the bounding box) as linear inequalities, supplying some linear objective function and passing everything to a linear programming (LP) solver. The only doubt is a large number of points on the grid and a worst case complexity of LP solver. So I am still looking for an algorithm of lower complexity or an explicit formula for such continuation.

Further information
I have a subroutine (in a large project on fluid dynamics) which accepts $C^0$ continuous functions $f(x,y)$ defined on $R^2$ and monotonous w.r.t. each argument. I have such functions defined only in a piece of $R^2$. I'm going to continue them to the whole $R^2$ as monotonous functions and pass them to the subroutine, to satisfy its requirements. I understand that there are many possible continuations, any of them will be suitable for me. What should happen in the worst case: the subroutine will return a solution outside the domain and I will use it as an indicator that solution inside the domain does not exist. Otherwise, if the function will not be continued to $R^2$, the subroutine will just reject the function and exit.

If the domain is the square $[x_0,x_1]\times[y_0,y_1]$, the following formula will do the job: $$f(x,y)=f(\min(\max(x,x_0),x_1),\min(\max(y,y_0),y_1))\\ +\left\{\begin{array}{l}x-x_0,\ x<x_0\\0, \ x_0\le x\le x_1\\x-x_1,x>x_1\end{array}\right.+\left\{\begin{array}{l}y-y_0,\ y<y_0\\0, \ y_0\le y\le y_1\\y-y_1,y>y_1\end{array}\right.$$

Concerning to LP-algorithm, the inequalities are $$f(x^{k+1},y^n)-f(x^k,y^n)\ge0,\ f(x^k,y^{n+1})-f(x^k,y^n)\ge0$$ for function values in points of the grid, while for the linear objective one can select the sum of l.h.s. of the inequalities in attempt to keep them as small as possible.

For strict monotony the above inequalities can be enforced to $$f(x^{k+1},y^n)-f(x^k,y^n)\ge\epsilon(x^{k+1}-x^k),\ f(x^k,y^{n+1})-f(x^k,y^n)\ge\epsilon(y^{n+1}-y^n)$$ where small positive constant $\epsilon$ sets a lower bound on derivatives.

As input, the subroutine takes a function pointer, so that we really can substitute continuous functions there. On the other hand, I have the function values sampled on a 2D grid and use bilinear interpolation in between.

The main question is how to continue (extrapolate) the function outside of its domain, preserving its monotony. Presumably at first on the grid till the bounding box of the domain with an unknown algorithm X, then till infinities using the formula above.

The formula does not represent a linear fit. It describes how to extrapolate a continuous function given on a square $[x_0,x_1]\times[y_0,y_1]$ to the whole plane as a function increasing w.r.t. each argument. It is obviously continuous, being represented as a sum of 3 continuous functions. $min$, $max$ are all continuous functions. Further terms are continuously joined at $x=x_{0,1}$, $y=y_{0,1}$. Similarly, it is monotonously increasing (w.r.t. each argument), being represented as a sum of 3 monotonously increasing functions. It is strictly increasing, since at least one function in the sum is strictly increasing, at any given point. It is not smooth, but here the smoothness is not required.