如果你的意思是质量$\int f(x,y)\; dxdy$在每个单元格上，那么加速它的最佳方法是并行化它。每个单元格都是独立的，所以你应该把外循环分成$p$零件并同时完成$p$处理器。你不能做得比这更好。你想要的需要 $O(h^{-2})$集成。

如果您不需要计算质量的单元格的统一网格，您可以研究自适应网格细化。如果您的地区有$f$如果事情非常好，您可能会通过不做不必要的工作而获得显着的加速。您这样做的代价是计算和使用结果的大量开销。$400\times400$可能太小而看不到自适应方法的好处。

最后，您在每个单元内使用什么集成方法？如果该方法是正交的，并且正交的点太多，您可能会受益于将其在单元格内的顺序减少到更合适的顺序。这只会降低复杂性的常数，但它可能是一个数量级或更多的差异。

在 Matlab 中嵌套循环可能很慢，通常消除循环的“矢量化”解决方案会更快。慢速循环曾经是 Matlab 中的常态，但自从 MathWorks 引入了他们的“即时编译器”以来，这种情况就不那么正确了。

但是，对于某些计算，矢量化解决方案仍然可以更快，因为许多内置的 Matlab 函数是多核的（即自动并行运行，当这样做很有效时）。

在不知道您的问题的确切细节的情况下，我不能肯定地说，但这是一种可能会有所帮助的方法。这本质上是有效计算移动平均值的“积分图像”方法，运行时间与滤波器大小无关。

首先，在一维空间中，假设你有一个网格$\mathbf{x}=[x_1,...x_{n+1}]$，以及矢量化函数，并且您想要计算在单元格上的积分，即。然后一个简单的方法是在单元边缘，并假设是分段线性的积分：$f=@(x)...$$f$$I_k=\int_{x_k}^{x_{k+1}}f(x)dx$$f$$\mathbf{f}=f(\mathbf{x})$$\mathbf{f}$

f=fnc(x); dx=diff(x); fmid=conv(f,[1,1],'valid')/2; I=fmid.*dx;

这只是梯形规则。

为了提高准确性，您可能需要对进行采样，而不仅仅是每个单元的端点。然后诀窍是进行采样，但对累积质量 (CDF) 进行子采样。这和上面一样，但略有变化：$f$$f$

F=cumsum([0,fmid.*dx]); I=diff(F(1:lag:end));

例如，“滞后”可能是 10。

在 2D 中，这种方法是这样的：

f=bsxfun(fnc,x,y'); dA=diff(x)*diff(y)'; fmid=conv2(f,ones(2),'valid')/4;
F=cumsum(cumsum(padarray(fmid.*dA,[1,1],'pre')),2);
I=diff(diff(F(1:lag:end,1:lag:end)),[],2);