一般情况下,对于任意一组多项式不等式,除了枚举搜索之外,没有其他方法可以解决这个问题,即使有有限数量的解也是如此。(或许我应该说不是不可能,因为枚举搜索可以解决这个问题,但需要注意浮点数问题。) 需要注意的是,对于高阶不等式,感兴趣的定义域不一定是单连通的。
编辑:原作者询问如何进行搜索。
考虑以下问题:
x^3 + y^3 >= 1e12
x^4 + y^4 <= 1e16
x >= 0, y >= 0
解决此系统的所有整数解。请注意,任何形式的整数规划都无法满足此处请求的所有整数解。
在这里使用网格会强制我们查看域(0:10000)X(0:10000)中的点。因此,它将强制我们对一组1e8个点进行采样,并测试每个点是否满足约束条件。
简单的循环可能比那更有效,尽管仍需要一些努力。
tic
xmax = 10000;
xy = cell(1,xmax);
for x = 0:xmax
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
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位机器上运行了这个测试,拥有8GB的内存。但即便如此,这个测试本身也是一个CPU占用量很高的程序。
Elapsed time is 50.182385 seconds.
请注意,浮点数计算有时会导致找到不同数量的点,具体取决于计算方法。
最后,如果您的约束方程更加复杂,您可能需要在y的边界表达式中使用根来帮助确定约束满足的位置。好处在于它仍然适用于更复杂的多项式边界。