Matlab 中多边形的双重积分

Question

我有一个函数@f(x,y)，我想在 MATLAB 中计算该函数在某个凸多边形上的积分。多边形不一定是矩形，这就是为什么我不能使用 MATLAB 的函数“dblquad”。我拥有的多边形由向量 X 和 Y 表示的一组顶点给出，即顶点为 (X(1),Y(1)),....,(X(n),Y(n) ）。有什么功能或方法可以使用吗？

Answer 1

诀窍是使用工具在感兴趣的区域内进行集成。我编写了一些用于在三角域中集成的工具。

% Define a function to integrate.
% This function takes an nx2 array, where each row
% contains a single point to evaluate the kernel at.
% This computes x^2 + y^2 at each point.
fun = @(xy) sum(xy.^2,2);

% define the domain as a triangulated polygon
% this tool uses ear clipping to do so.
sc = poly2tri([1 4 3 1],[1 3 5 4]);

% Gauss-Legendre integration over the 2-d domain
[integ,fev]= quadgsc(fun,sc,2)
integ =
       113.166666666667
fev =
     8

% the triangulated polygon...
plotsc(sc,'facecolor','none','markerfacecolor','r')
axis equal
grid on

我们可以将函数本身可视化为该多边形域上的映射 z(x,y)。当提供范围字段时，单纯复形将变成来自 2-d (x,y) 域的 2-1 映射。

sc2 = refinesc(sc,'max',.5);
sc2.range = fun(sc2.domain);
plotsc(sc2,'markerfacecolor','r')
grid on
view(17,12)

这是感兴趣域上的简单多项式函数，因此默认的低阶高斯积分就足够了。使用的方案是三角形张量积形式的高斯-勒让德方案，不是真正最优的，但可行。高斯求积的问题是它不具有适应性。它基于有限点集上多项式的隐式近似来计算估计值。

上述估计使用 8 个函数评估来计算该估计。由于内核是一个低阶多项式，因此它应该表现得很好。问题是，您需要知道这是否是正确的解决方案。这是高斯求积的问题，没有简单的方法可以知道答案是否正确，除非用更高阶的方案解决问题，直到它看起来收敛。

看到，重心处每个三角形有 1 个点，我们得到了错误的答案，但高阶估计都一致。

[integ,fev]= quadgsc(fun,sc,1)
integ =
       107.777777777778
fev =
     2

[integ,fev]= quadgsc(fun,sc,3)
integ =
       113.166666666667

fev =
    18

[integ,fev]= quadgsc(fun,sc,4)
integ =
       113.166666666667
fev =
    32

编写完quadgsc后，我不得不尝试一个自适应求解器，它的工作方式与MATLAB中其他四边形工具的工作方式相同。这会对三角测量进行自适应细化，寻找解不稳定的三角形。问题是，我从来没有以令我满意的方式完成这些工具的编写。对于三角域上的体积问题，可以采用许多不同的方法。 Quadrsc 进行低阶解，然后对其进行细化，使用理查森外推法，然后比较结果。对于差异太大的任何三角形，它会进一步细化它们，直到收敛。

例如，

[integ,fev]= quadrsc(fun,sc)
integ =
          113.166666666667

fev =
    16

所以这有效。这个问题出现在更复杂的内核上，其中的问题是知道何时停止细化，并在用完太多函数评估之前停止细化。我从来没有完全让我满意，所以我从来没有发布过这些工具。我可以将工具箱发送给那些直接向我发送邮件的人。该 zip 文件大约 2.4 MB。有一天我会抽出时间来完成这些工具，我希望......