Finden Sie die Discrete Paar {x, y}, dass Satisfy Inequality Constriants

Question

Ich habe einige Ungleichheiten in Bezug auf {x,y}, daß die folgenden Gleichungen erfüllt:

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

Beachten Sie, dass x und y integer sein müssen.

Grafisch kann es wie folgt dargestellt werden, der blaue Bereich ist der Bereich, der die obigen Ungleichungen erfüllt:

Die Frage ist nun, gibt es eine Funktion in Matlab, die jede zulässige Paar findet von {x,y}? Wenn es einen Algorithmus gibt, um so etwas zu tun, würde ich mich freuen, davon zu hören.

Natürlich können wir immer einen Brute-Force-Ansatz verwenden, bei dem wir jede mögliche Kombination von {x,y} testen, um festzustellen, ob die Ungleichungen erfüllt sind. Aber das ist der letzte Ausweg, weil es zeitaufwendig ist. Ich suche nach einem cleveren Algorithmus, der dies tut, oder im besten Fall, einer vorhandenen Bibliothek, die ich sofort verwenden kann.

Die x^2+y^2>=100and x^2+y^2<=200 sind nur Beispiele; in Wirklichkeit können f und g beliebige Polynomfunktionen eines beliebigen Grades sein.

Bearbeiten: C# -Code sind auch willkommen.

Answer 1

Dies ist sicherlich nicht möglich, für einen allgemeinen Satz von Polynom Ungleichheiten im Allgemeinen zu tun, durch ein anderes Verfahren als enumerative suchen, auch wenn es eine endliche Anzahl von Lösungen. (Vielleicht sollte ich sagen, nicht trivial, wie es möglich ist. Enumerative Suche funktioniert, abhängig von Gleitkomma-Problemen.) Beachten Sie, dass der Interessenbereich nicht einfach für Ungleichungen höherer Ordnung verbunden sein muss.

Bearbeiten: Das OP hat gefragt, wie man eine Suche durchführen könnte.

Betrachten Sie das Problem

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

x >= 0, y >= 0 

für alle ganzzahligen Lösungen dieses Systems lösen. Beachten Sie, dass Integer-Programmierung in JEDER Form hier nicht ausreicht, da ALLE Integer-Lösungen angefordert werden.

Die Verwendung von Meshgrid würde uns zwingen, Punkte in der Domäne (0: 10000) X (0: 10000) zu betrachten. Es würde uns also zwingen, eine Reihe von 1e8-Punkten zu testen und jeden Punkt zu testen, um zu sehen, ob sie die Beschränkungen erfüllen.

Eine einfache Schleife kann möglicherweise effizienter sein, obwohl es immer noch einige Anstrengungen erfordern wird.

% 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 

Die Zeit, die erforderlich war ...

Elapsed time is 0.600419 seconds. 

Von den 100.020.001 Kombinationen, die wir für getestet haben könnte, wie viele Lösungen fanden wir?

size(xy) 
ans = 
      2  4371264 

Zugegeben, die erschöpfende Suche ist einfacher zu schreiben.

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

I lief dies auf einer 64-Bit-Maschine mit 8 gig RAM. Aber selbst der Test selbst war ein CPU-Schwein.

Elapsed time is 50.182385 seconds. 

Beachten Sie, dass Gleitkommaüberlegungen manchmal dazu führen, dass je nach Ausführung der Berechnungen eine andere Anzahl von Punkten gefunden wird.

Wenn Ihre Abhängigkeitsbedingungen komplexer sind, müssen Sie möglicherweise im Ausdruck für die Grenzen von y Wurzeln verwenden, um festzustellen, wo die Bedingungen erfüllt sind. Das Schöne daran ist, dass es immer noch für kompliziertere Polynomgrenzen funktioniert.