如何在MATLAB中对离散二维曲面进行积分集成？

Question

我有一个函数 z = f(x, y)，其中 z 是点 (x, y) 的值。如何在 MATLAB 中对 x-y 平面上的 z 进行积分？

上述函数实际上类似于哈希表。也就是说，给定一个 (x, y) 对，我可以查找表格以找到相应的 z 值。

如果点在 x-y 平面上均匀分布，则问题会非常简单，此时可以将所有 z 值相加，乘以底部面积，并最终除以点数即可。但是，分布不均匀，如下所示。因此，我实际上正在寻求最小化误差的计算方法。

- Sibbs Gambling

---

那么你知道函数定义的所有 x-y 值吗？ - shimizu

@shimizu 是的。基本上，对于每一组x和y，我都有一个相应的z。我知道它们全部。上面的图是（x，y）的散点图。 - Sibbs Gambling

那么你的问题本质上是关于在非均匀的 x 和 y 集合中进行 2D 积分吗？ - shimizu

@shimizu 是的，没错。 - Sibbs Gambling

Answer 1

目前被接受的答案只适用于网格数据。如果您的数据是分散的，您可以使用以下方法：

scatteredInterpolant + integral2：

f = scatteredInterpolant(x(:), y(:), z(:), 'linear');
int = integral2(@(x,y) f(x,y), xmin, xmax, ymin, ymax);

这定义了数据z(i) = f(x(i),y(i))的线性插值函数f，并将其用作integral2的参数。请注意，ymin和ymax可以是依赖于x的函数句柄，而不是双精度浮点数。因此，通常您将集成矩形，但这可以用于集成稍微复杂一些的区域。

如果您的积分区域比较复杂或具有空洞，则应考虑对您的数据进行三角剖分。

自己动手使用triangulation：

假设您的积分区域由三角剖分trep给出，例如可以通过trep = delaunayTriangulation(x(:), y(:))获得。如果您的值z与z(i) = f(trep.Points(i,1), trep.Points(i,2))相对应，您可以使用以下集成例程。它计算线性插值的精确积分。这是通过计算所有三角形的面积，然后将这些面积用作每个三角形上的中点（平均）值的权重来完成的。

function int = integrateTriangulation(trep, z)
P = trep.Points; T = trep.ConnectivityList;
d21 = P(T(:,2),:)-P(T(:,1),:);
d31 = P(T(:,3),:)-P(T(:,1),:);
areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
int = areas'*mean(z(T),2);

- knedlsepp

---

感谢三角测量程序。如果您有一个非常不均匀的数据集（例如具有显着角落细化的矩形区域），则此程序非常有用。如果您强制该细化在整个网格上传播，则强制实施规则矩形网格的方法可能会变得非常低效和内存密集。 - Nick J

Answer 2

如果你有一个离散数据集，其中你拥有定义z的所有x和y值，那么只需获取与这些（x，y）对应的Zdata矩阵。保存此矩阵，然后您可以使用interp2将其变成连续函数：

function z_interp = fun(x,y)

    z_interp = interp2(Xdata,Ydata,Zdata,x,y);

end

然后你可以使用integral2函数来计算积分：

q = integral2(@fun,xmin,xmax,ymin,ymax)

其中@fun是一个函数句柄，它接受两个输入参数。

- shimizu

---

问题是，我没有一个函数。相反，我有一系列离散的点。 - Sibbs Gambling

那么，你有一个定义了 z 的 x 和 y 点的列表吗？它们是均匀分布的吗？ - shimizu

非常感谢，请您查看更新后的问题。 - Sibbs Gambling

就快完成了，除了我遇到了“积分被积函数输出大小与输入大小不匹配”的错误。也许这是因为我的积分被积函数的输出是一个复数？ - Sibbs Gambling

Answer 3

最近我在MatLab中需要集成一个双变量正态分布。这个想法非常简单。Matlab通过网格定义了一个表面，因此从x、y开始，你需要这样做：

x = -10:0.05:10;
y = x;

[X,Y] = meshgrid(x',y');

...例如。然后，让我们称FX为定义表面每个点值的函数。要计算积分，您只需要执行以下操作：

surfint = zeros(length(X),1);
for a = 1:length(X)
   surfint(a,1) = trapz(x,FX(:,a));     
end

trapz(x, surfint)

对我来说，这是最简单的方法。

- Paulo Roberto C Passos