trova la discreta Coppia di {x, y} che soddisfano disuguaglianza Constriants

Question

ho un paio di disuguaglianze in materia di {x,y}, che soddisfa le seguenti equazioni:

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

noti che x e y deve essere intero. ci

La domanda è, è qualsiasi funzione in Matlab che trova ogni coppia ammissibile:

Graficamente può essere rappresentata come segue, la regione blu è la regione che soddisfa le disuguaglianze sopra di {x,y}? Se c'è un algoritmo per fare questo genere di cose sarei lieto di sentirlo anche a questo proposito.

Ovviamente, un approccio che possiamo sempre utilizzare è l'approccio a forza bruta in cui testiamo ogni possibile combinazione di {x,y} per verificare se le disuguaglianze sono soddisfatte. Ma questa è l'ultima risorsa, perché richiede molto tempo. Sto cercando un algoritmo intelligente che faccia questo, o nel migliore dei casi, una libreria esistente che posso usare direttamente.

x^2+y^2>=100and x^2+y^2<=200 sono solo esempi; in realtà f e g possono essere funzioni polinomiali di qualsiasi grado.

Modifica: codice C# sono i benvenuti.

Answer 1

Ciò non è sicuramente possibile fare in generale per un insieme generale di disuguaglianze polinomiali, con qualsiasi metodo diverso dalla ricerca enumerativa, anche se esiste un numero finito di soluzioni. (Forse dovrei dire non banale, come è possibile. La ricerca enumerativa funzionerà, soggetta a problemi in virgola mobile.) Si noti che il dominio di interesse non deve essere semplicemente connesso per disuguaglianze di ordine superiore.

Modifica: l'OP ha chiesto come procedere per effettuare una ricerca.

consideri il problema

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

x >= 0, y >= 0 

risolvere per tutte le soluzioni intere di questo sistema. Nota che la programmazione intera in QUALSIASI forma non è sufficiente qui, poiché sono richieste TUTTE le soluzioni a numeri interi.

L'uso di meshgrid qui ci costringerebbe a guardare i punti nel dominio (0: 10000) X (0: 10000). Quindi ci costringerebbe a campionare un set di 1e8 punti, testando ogni punto per vedere se soddisfano i vincoli.

Un ciclo semplice può potenzialmente essere più efficiente di quello, anche se richiederà ancora un certo sforzo.

% 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 

Il tempo necessario è stato ...

Elapsed time is 0.600419 seconds. 

Delle 100020001 combinazioni che potremmo avere testati per, quante soluzioni abbiamo trovato?

size(xy) 
ans = 
      2  4371264 

Certamente, la ricerca esauriente è più semplice da scrivere.

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

L'ho eseguito su una macchina a 64 bit, con 8 gig di ram. Ma anche così il test stesso era un maiale da CPU.

Elapsed time is 50.182385 seconds. 

nota che le considerazioni virgola mobile a volte causare un diverso numero di punti da trovare, a seconda di come i calcoli verranno eseguiti.

Infine, se le equazioni dei vincoli sono più complesse, potrebbe essere necessario utilizzare le radici nell'espressione per i limiti su y, per identificare dove sono soddisfatti i vincoli. La cosa bella qui è che funziona ancora per limiti polinomiali più complicati.