Znajdź dyskretne Para {x, y}, które spełniają nierówność Constriants

Question

Mam kilka nierówności dotyczących {x,y}, który spełnia następujące równania:

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

Zauważ, że x i y musi być liczbą całkowitą. nie

Teraz powstaje pytanie, jest dowolną funkcją w Matlab, która znajduje każdą dopuszczalną pary:

Graficzne mogą być przedstawione jak następuje region niebieski jest regionem, który spełnia powyższe nierówności z {x,y}? Jeśli istnieje algorytm do robienia tego rodzaju rzeczy, chętnie bym się o tym dowiedział.

Oczywiście jednym podejściem, które zawsze możemy zastosować, jest podejście brutalnej siły, w którym testujemy każdą możliwą kombinację {x,y}, aby sprawdzić, czy nierówności są spełnione. Ale to jest ostatnia deska ratunku, bo to czasochłonne. Szukam sprytnego algorytmu, który to robi, lub w najlepszym razie istniejącej biblioteki, której mogę użyć od razu.

The x^2+y^2>=100and x^2+y^2<=200 to tylko przykłady; w rzeczywistości f i g może być dowolną funkcją wielomianową dowolnego stopnia.

Edytuj: Kod C# również jest mile widziany.

Answer 1

Zasadniczo nie da się zrobić ogólnego zestawienia nierówności wielomianowych, za pomocą dowolnej metody, innej niż wyszukiwanie enumeratywne, nawet jeśli istnieje skończona liczba rozwiązań. (Być może powinienem powiedzieć, że nie jest to banalne, ponieważ jest to możliwe, ponieważ wyszukiwanie wyliczeniowe będzie działało, z zastrzeżeniem kwestii zmiennoprzecinkowych.) Należy zauważyć, że domena zainteresowania nie musi być po prostu połączona w przypadku nierówności wyższego rzędu.

Edytuj: OP zadał pytanie, w jaki sposób można przystąpić do wyszukiwania.

rozważyć problem

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

x >= 0, y >= 0 

rozwiązać dla wszystkich rozwiązań całkowitych tego systemu. Zauważ, że całkowite programowanie w DOWOLNEJ formie nie wystarczy tutaj, ponieważ żądane są WSZYSTKIE rozwiązania liczb całkowitych.

Użycie siatki siatkowej zmusiłoby nas do spojrzenia na punkty w domenie (0: 10000) X (0: 10000). Zmusiłoby to nas do wypróbowania zestawu 1e8 punktów, testowania każdego punktu, aby sprawdzić, czy spełniają one ograniczenia.

Prosta pętla może być potencjalnie bardziej wydajna, chociaż nadal będzie wymagać pewnego wysiłku.

% 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 

Czas potrzebny był ...

Elapsed time is 0.600419 seconds. 

Spośród 100020001 kombinacji, które moglibyśmy testowanych, ile rozwiązań nie znajdziemy?

size(xy) 
ans = 
      2  4371264 

Oczywiście, wyczerpujące wyszukiwanie jest łatwiejsze do napisania.

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

Uruchomiłem to na komputerze 64-bitowym, z 8 gig pamięci RAM. Ale nawet sam test był wieżą CPU.

Elapsed time is 50.182385 seconds. 

Należy pamiętać, że rozważania zmiennoprzecinkowe czasami powodują różną liczbę punktów, w zależności od sposobu obliczeń.

Na koniec, jeśli równania ograniczenia są bardziej złożone, może być konieczne użycie pierwiastków w wyrażeniu dla granic y, aby pomóc określić, gdzie są spełnione ograniczenia. Przyjemnie jest, że wciąż działa na bardziej skomplikowane wielomianowe granice.