Encontrar el par de discretos {x, y} que satisfacen la desigualdad Restricciones

Question

Tengo algunas desigualdades en relación con {x,y}, que satisface las siguientes ecuaciones:

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

Nota que x y y deben estar entero.

Gráficamente se puede representar de la siguiente manera, la región azul es la región que satisface las desigualdades anteriores:

La pregunta ahora es, ¿hay alguna función en Matlab que encuentra cada par admisible de {x,y}? Si hay un algoritmo para hacer este tipo de cosas que se alegraría de oír hablar de eso también.

Por supuesto, uno de los enfoques que podamos siempre uso es método de fuerza bruta, donde ponemos a prueba todas las combinaciones posibles de {x,y} para ver si las desigualdades son satisfechas. Pero este es el último recurso, debido a consumirlo de tiempo. Estoy buscando un algoritmo inteligente que hace esto, o en el mejor de los casos, una biblioteca existente que puedo usar recta de distancia.

El x^2+y^2>=100 and x^2+y^2<=200 son ejemplos sólo; en realidad f y g puede ser cualquiera de las funciones polinómicas de cualquier grado.

Editar:. C # código son bienvenidos, así

Answer 1

Esto no es, sin duda posible hacerlo, en general, para un conjunto general de las desigualdades polinomiales, por cualquier método que no sea la búsqueda enumerativa, incluso si hay un número finito de soluciones. (Tal vez debería decir no es trivial, ya que es posible. Enumerativo de búsqueda va a funcionar, sujetos a cuestiones de punto flotante.) Tenga en cuenta que el dominio de la necesidad interés no puede conectar simplemente por las desigualdades de orden superior.

Editar:. El PO ha preguntado sobre cómo se puede proceder a hacer una búsqueda

Considere el problema

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

x >= 0, y >= 0

Resolver para todas las soluciones enteros de este sistema. Tenga en cuenta que la programación entera en ninguna forma no será suficiente aquí, ya que todos los enteros se solicitan soluciones.

El uso de meshgrid aquí nos obligaría a mirar puntos en el dominio (0: 10000) X (0: 10.000). Por lo que nos obligaría a la muestra un conjunto de puntos 1e8, poniendo a prueba todos los puntos para ver si cumplen las restricciones.

Un simple bucle potencialmente puede ser más eficiente que, aunque todavía se requiere un cierto esfuerzo.

% 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

El tiempo requerido era ...

Elapsed time is 0.600419 seconds.

De los 100020001 combinaciones que podríamos haber probado para, cuántas soluciones Qué encontramos?

size(xy)
ans =
           2     4371264

Es cierto que la búsqueda exhaustiva es más fácil de escribir.

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

me encontré con este en una máquina de 64 bits, con 8 GB de RAM. Pero aún así la prueba en sí era un cerdo de la CPU.

Elapsed time is 50.182385 seconds.

Tenga en cuenta que las consideraciones de punto flotante veces causarán un número diferente de puntos que se encuentran, dependiendo de cómo se hacen los cálculos.

Por último, si sus ecuaciones de restricción son más complejos, puede que tenga que utilizar raíces en la expresión de los límites de y, para ayudar a identificar cuando se cumplan las restricciones. Lo bueno aquí es que aún trabaja para los límites polinómicas más complicado.