Найдите дискретные пары {x, y}, которые удовлетворяют неравенству

Question

У меня есть несколько неравенств в отношении {x,y}, что удовлетворяет следующим уравнениям:

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

Обратите внимание, что x а также y должен быть целым числом.

Графически его можно представить следующим образом, Голубой регион является областью, которая удовлетворяет вышеуказанному неравенством:

Вопрос сейчас есть, есть ли функция в Matlab, которая находит каждую допустимую пару {x,y}? Если есть алгоритм, чтобы сделать так, я был бы рад слышать об этом.

Конечно, один подход, который мы всегда можем использовать, - это грузоподъемность, в котором мы тестируем все возможные комбинации {x,y} чтобы увидеть, удовлетворены ли неравенство. Но это последний курорт, потому что поглощение времени. Я ищу умного алгоритм, который делает это или в лучшем случае, существующая библиотека, которую я могу использовать прямо.

То x^2+y^2>=100 and x^2+y^2<=200 просто примеры; в реальности f а также g Может быть любые полиномиальные функции любой степени.

Редактировать: C # код также приветствуется.

Answer 1

Это, безусловно, невозможно сделать в целом для общего набора полиномиальных неравенств, по любому способу, отличным от перечисленного поиска, даже если есть конечное количество решений. (Возможно, я должен сказать, что не тривиально, как это возможно. Перечислительный поиск будет работать, при условии, что проблемы с плавающей точкой.) Обратите внимание, что домен интересов не должен быть просто подключен к неравенству более высокого порядка.

Редактировать: OP задал вопрос о том, как можно приступить к поиску.

Рассмотрим проблему

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

x >= 0, y >= 0

Решить для всех целочисленных решений этой системы. Обратите внимание, что целочисленное программирование в любой форме здесь не хватит, поскольку все целочисленные решения запрашиваются.

Использование MeshGrid здесь заставило бы нас взглянуть на точки домена (0: 10000) x (0: 10000). Таким образом, это заставило бы нас попробовать набор очков 1E8, тестируя каждую точку, чтобы увидеть, удовлетворяют ли они ограничениям.

Простая петля может быть более эффективной, чем это, хотя все равно потребует некоторых усилий.

% 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

Требуется время ...

Elapsed time is 0.600419 seconds.

Из 100020001 комбинаций, на которые мы могли бы проверить, сколько решений мы нашли?

size(xy)
ans =
           2     4371264

По общему признанию, исчерпывающий поиск проще для записи.

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

Я провел это на 64-битной машине, с 8 гигом оперативной памяти. Но даже так сами тест был боров CPU.

Elapsed time is 50.182385 seconds.

Обратите внимание, что соображения с плавающей точкой иногда вызывают другое количество точек, которые необходимо найти, в зависимости от того, как выполняются вычисления.

Наконец, если ваши уравнения ограничения более сложны, вам может потребоваться использовать корни в выражении для границ на Y, чтобы помочь определить, где ограничения удовлетворяются. Приятно здесь это все еще работает для более сложных полиномиальных границ.