Find a discrete pair {x, y} that satisfies the inequalities

Question

I have several inequalities with respect to {x,y} that satisfy the following equations:

 x>=0 y>=0 f(x,y)=x^2+y^2>=100 g(x,y)=x^2+y^2<=200 

Note that x and y must be integers.

Graphically, this can be represented as follows: the blue region is the region satisfying the above inequalities:

The question is, is there any function in Matlab that finds every valid pair {x,y} ? If there is such an algorithm, I would also be glad to hear about it.

Of course, one approach we can always use is brute force, where we check every possible combination of {x,y} to see if the inequalities are true. But this is the last resort, because it takes a lot of time. I am looking for a smart algorithm that does this or, at best, an existing library that I can use right now.

x^2+y^2>=100 and x^2+y^2<=200 - only examples; in fact, f and g can be any polynomial functions of any degree.

Edit: C # code is also welcome.

Answer 1

Of course, this is generally impossible to do for a general set of polynomial inequalities using any method other than an enumerable search, even if there are a finite number of solutions. (Perhaps I should say that this is not trivial as possible). An enumeration search will work with floating point problems.) Note that the area of ​​interest does not have to be simply related to higher-order inequalities.

Edit: The OP asked a question about how to continue the search.

Consider the problem

 x^3 + y^3 >= 1e12 x^4 + y^4 <= 1e16 x >= 0, y >= 0 

Solve for all integer solutions of this system. Please note that integer programming in ANY form is not enough here, since ALL whole solutions are required.

Using meshgrid here will make us look at the points in the domain (0: 10000) X (0: 10000). Thus, this will force us to try a set of 1e8 points, checking each point to see if they satisfy the constraints.

A simple cycle could potentially be more efficient than this, although it would still require some effort.

 % Note that I will store these in a cell array, % since I cannot preallocate the results. tic xmax = 10000; xy = cell(1,xmax); for x = 0:xmax % solve for y, given x. This requires us to % solve for those values of y such that % y^3 >= 1e12 - x.^3 % y^4 <= 1e16 - x.^4 % These are simple expressions to solve for. y = ceil((1e12 - x.^3).^(1/3)):floor((1e16 - x.^4).^0.25); n = numel(y); if n > 0 xy{x+1} = [repmat(x,1,n);y]; end end % flatten the cell array xy = cell2mat(xy); toc 

The required time was ...

 Elapsed time is 0.600419 seconds. 

Of the 100,000,000 combinations that we could test, how many solutions did we find?

 size(xy) ans = 2 4371264 

Admittedly, an exhaustive search is easier to write down.

 tic [x,y] = meshgrid(0:10000); k = (x.^3 + y.^3 >= 1e12) & (x.^4 + y.^4 <= 1e16); xy = [x(k),y(k)]; toc 

I ran this on a 64 bit machine with 8 gigabytes. But even the test itself was a CPU failure.

 Elapsed time is 50.182385 seconds. 

Note that floating point considerations sometimes cause a different number of points, depending on how the calculations are performed.

Finally, if your constraint equations are more complex, you may need to use the roots in the expression for the boundaries on y to determine where the constraints are satisfied. It's nice that this still works for more complex polynomial boundaries.