Найдите дискретную пару из {x, y}, удовлетворяющих константам неравенства - PullRequest
Новая
       43

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

5 голосов
/ 25 сентября 2010

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

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

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

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

alt text

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

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

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

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

1 Ответ

4 голосов
/ 25 сентября 2010

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

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

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

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

x >= 0, y >= 0

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

Использование сетки здесь заставит нас посмотреть на точки в домене (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 гигабайтами оперативной памяти. Но даже в этом случае сам тест был загружен процессором. 

Elapsed time is 50.182385 seconds.

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

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

...