紧支撑函数在三角形上的数值积分

计算科学 正交
2021-12-19 06:10:06

正如标题所示,我正在尝试计算三角形上紧支持函数(温德兰五次多项式)的积分。请注意,函数的中心位于 3-D 空间中的某个位置。我将此函数集成在一个任意但很小的三角形()上。我目前正在使用 Dunavant,1985 (p=19) 描述的集成。area<(radius/4)22

然而,这些求积规则似乎不适用于紧支持的问题。这得到了以下事实的支持:当我在使用三角形离散化的平面上积分(半径为 1 的圆内为 1 的函数)时,我的(归一化)结果介于1.001 和 0.897。f(r)=[r1]

所以我的问题是,这种问题是否存在专门的求积法则?低阶复合积分规则会更好吗?

不幸的是,这个例程在我的代码中非常关键,因此精度至关重要。另一方面,我需要为单个时间步执行“几次”这种集成,因此计算成本不应该太高。并行化不是问题,因为我将串行执行集成本身。

提前感谢您的回答。

编辑:Wendland 的五次多项式由其中r_0是\mathbb{R}^3中的任意向量W(q)=[q2]αh3(1q2)4(2q+1)α=2116πq=rr0hr0R3

EDIT2:如果Δ是二维三角形,那么我想\omega(r) = W(\frac{\|r-r_0\|}{h})计算\int_\Delta \omega(r) dr . 所以W中的q永远不会小于 0。请注意,积分是在\mathbb{R}^3中的二维曲面上的曲面积分Δω(r)drω(r)=W(rr0h)qWR3

EDIT3:我有一个一维(线)问题的分析解决方案。也可以为二维(三角形)计算一个。

3个回答

由于函数在内是平滑的,但不是固定度数(即在平面内),我建议在两个维度上使用简单的自适应方案,例如带有Romberg 方法的梯形规则。q2

也就是说,如果您的三角形由顶点定义,并且您有一个沿线从to积分的例程,您可以执行以下操作(在 Matlab 表示法中):xyzR3romb(f,a,b)fab

int = romb( @(xi) romb( W , xi , y+(z-y)*(xi-x)./(z-x) ) , x , z );

romb中,不要使用固定数量的点,而是不断扩大表格,直到两个连续对角线之间的差值低于您所需的容差。由于您的函数是平滑的,因此这应该是一个很好的误差估计。

的域之外,您可以尝试相应地调整上述代码中的积分限制。W(q)

这可能不是解决问题的计算效率最高的方法,但与固定度数规则相比,自适应性会给您带来更多的鲁棒性。

有关容积规则的详细概述,请参阅“R. Cools, An Encyclopaedia of Cubature Formulas J. Complexity, 19: 445-453, 2003”。使用固定规则可以为您提供一些规则精确积分多项式的优势(就像高斯求积在一维中所做的那样)。

Cools 也是 CUBPACK 的主要作者之一,CUBPACK是一个用于数值容积的软件包。

积分规则假设该函数在局部被一个低次多项式很好地逼近。您的问题与紧凑型支持无关。紧支撑的径向基函数在支撑边界处是光滑的,并且可以毫无问题地使用高达光滑量级的正交规则。(高阶规则没有帮助;因此您可能不应该使用精确积分 5 次多项式的规则。)

在您的情况下,不准确性来自这样一个事实,即在您的情况下,对于附近的三角形,良好多项式近似性的假设失败,即使它们不包含r0r0

W的函数是平滑的,但中变为无限大积分超过,复合函数是的非光滑函数。qqrrr0rr

如果三角形不包含,则函数为但这无济于事,因为高阶导数在附近增长得非常快,并且高阶方法的误差与高阶导数成正比,因此非常大!r0Cinfr0

简单的补救方法是将每个三角形 T 分成 N_T 个子三角形。您可以将远离,并将靠近您可以离线计算出对于给定直径和距必须有多大才能达到所需的精度。的低阶公式NT=1r0NT1r0NTr0r0

当您在三角形上积分时,但是 3 维的,三角形显然在中。r0R3

的积分制表为三角坐标的函数(通过将其旋转到二维平面中进行归一化,使得一个顶点位于轴上,并将其反映为第二个顶点位于其上方)。该表格必须足够详细,以使线性或二次插值足够准确。但是您可以使用首先概述的慢速方法来创建此表。r0=0xyx

摆脱这个问题的另一种方法是使用紧凑支持的径向基函数,它是而不是中的多项式。这在任何地方都很流畅,并且易于集成。q2q