空间三角形上的数值积分

The Numerical Integration on Spatial Triangle

摘要: 本文通过引进体积坐标，结合数值积分，给出了统一求空间三角形上的数值积分方法。该方法方便快捷，易于程序化。最后本文给出了用该方法求任意函数数值积分的MATLAB程序。

Abstract: This paper presents a method to solve numerical integration on spatial triangle by combing volume coordinates and numerical integration. This method is easy and efficient and can be programed. Finally, this paper gives the MATLAB program for any function to solve the numerical integration.

1. 引言

在用有限元等数值方法计算求解工程中的数学问题时，我们通常会采用四面体(每个四面体我们称之为单元)对空间区域进行剖分。不可避免的，比如在计算残差项的时候，我们通常需要在四面体的一个面(这个面是一个空间三角形)上进行积分。理论上，我们可用第一型曲面积分公式 [1] ：

${\displaystyle {\iint }_{S}f\left(x,y,z\right)\text{d}S}={\displaystyle {\iint }_{D}f\left(x,y,z\left(x,y\right)\right)\sqrt{1+z_x^2+z_y^2}\text{d}x\text{d}y}$ (1)

进行计算。但是，首先，四面体的各个面的大小、形状不一，统一求积分必然是一件困难的事情。其次被积函数有可能比较复杂，求得积分公式难以求得准确解，迫使我们不得不转向数值积分。最后，若用数值积分，将Oxy平面上单位三角形上的积分节点转化为空间三角形上的积分节点也不是件容易的事。

本文提供的方法，很好的解决了上述三个问题。使得我们在实际应用中能够快捷方便的求解任意一个函数在任意一个空间三角形上的积分。

2. 求积方法的思想

2.1. 空间直角坐标与体积坐标

我们设一个空间三角形ABC的三个顶点在空间直角坐标系下的坐标为 $A\left(x_1,y_1,z_1\right)$ ， $B\left(x_2,y_2,z_2\right)$ ， $C\left(x_3,y_3,z_3\right)$ ，为了引进体积坐标，我们需要引入不在 $\Delta ABC$ 所在平面的第四个点 $D\left(x_4,y_4,z_4\right)$ ，进而引入面积坐标 [2] ：

$p_1\left(x,y,z\right)=\left[\begin{array}{cccc}x& y& z& 1\\ x_2& y_2& z_2& 1\\ x_3& y_3& z_3& 1\\ x_4& y_4& z_4& 1\end{array}\right]/V$ , $p_2\left(x,y,z\right)=\left[\begin{array}{cccc}{x}_1& y_1& z_1& 1\\ x& y& z& 1\\ x_3& y_3& z_3& 1\\ x_4& y_4& z_4& 1\end{array}\right]/V$ ,

$p_3\left(x,y,z\right)=\left[\begin{array}{cccc}{x}_1& y_1& z_1& 1\\ x_2& y_2& z_2& 1\\ x& y& z& 1\\ x_4& y_4& z_4& 1\end{array}\right]/V$ , $p_4\left(x,y,z\right)=\left[\begin{array}{cccc}{x}_1& y_1& z_1& 1\\ x_2& y_2& z_2& 1\\ x_3& y_3& z_3& 1\\ x& y& z& 1\end{array}\right]/V$ .

其中 $$ 为四面体ABCD的体积。

因此，对于空间直角坐标系中任意一点 $p$ 的坐标 $\left(x,y,z\right)$ ，可以用体积坐标表示如下：

$\begin{array}{l}x=x_1{p}_1+x_2{p}_2+x_3{p}_3+x_4{p}_4,\\ y=y_1{p}_1+y_2{p}_2+y_3{p}_3+y_4{p}_4,\\ z=z_1{p}_1+z_2{p}_2+z_3{p}_3+z_4{p}_4.\end{array}$ (2)

性质1 [2] ： $p_i\left(x_j,y_j,z_j\right)=\{\begin{array}{l}0,i\ne j\\ 1,i=j\end{array},\left(i,j=1,2,3,4\right)$ 。

性质2 [2] ： $p_1+p_2+p_3+p_4=1$ 。

2.2. 对函数求积分

设函数 $u=u\left(x,y,z\right)$ 在空间三角形ABC有定义。由性质1知， $p_4$ 在∆ABC所在平面恒为零，所以 $p_1+p_2+p_3=1$ 。所以Oxy平面上直角边为1的等腰直角三角形上的积分节点 $\left(p_1,p_2\right)$ 容易转化为其上的体积积分节点 $\left(p_1,p_2,1-p_1-p_2\right)$ ，进而通过公式(2)映射成一般单元上的积分节点。

设 $\left(p_{1i},p_{2i},p_{3i}\right),i=1,2,\cdots ,n$ 是一组积分节点，相应积分节点的权记为 $w_i$ ，于是由公式(2)我们可以求得一般空间三角形上的积分节点 $\left(x_{1i},y_{2i},z_{3i}\right),i=1,2,\cdots ,n$ 。于是

${\displaystyle {\iint }_{\Delta ABC}u\left(x,y,z\right)}\approx S_{\Delta ABC}{\displaystyle {\sum }_i^{n}f\left(x_{1i},y_{2i},z_{3i}\right)}$ (3)

其中 $S_{\Delta ABC}$ 是∆ABC的面积。

3. 算法

Step 1：首先选择Oxy平面上单位三角形上的一组积分节点 $\left(p_{1i},p_{2i}\right),i=1,2,\cdots ,n$ 和相应的权 $w_i$ ，计算相应体积坐标积分节点 $\left(p_{1i},p_{2i},1-p_{1i}-p_{2i}\right),i=1,2,\cdots ,n$ 。

Step 2：由公式(2)计算空间三角形上的积分节点 $\left(x_{1i},y_{2i},z_{3i}\right),i=1,2,\cdots ,n$ 。

Step 3：用公式(3)计算积分。

4. 示例

例1： ${\displaystyle {\iint }_{S}xyz\text{d}s}$ ， $S$ 为平面 $x+y+z=1$ 在第一卦限的部分。

解法一：用公式(1)求解

$z=1-x-y,z_x=-1,z_y=-1,\text{d}s=\sqrt{3}\text{d}x\text{d}y$

所以

${\displaystyle {\iint }_{S}xyz\text{d}s}={\displaystyle {\iint }_{D_{xy}}xy\left(1-x-y\right)\sqrt{3}\text{d}x\text{d}y}={\displaystyle {\int }_0^1\text{d}x}{\displaystyle {\int }_0^{1-x}\sqrt{3}xy\left(1-x-y\right)\text{d}y}=\frac{\sqrt{3}}{120}$ .

解法二：本文参考 [3] ，选用具有三次代数精度的求积节点，并通过MATLAB编程(见附录)求解如下：

程序输入：node=[1,0,0;0,1,0;0,0,1];

elem=[1,2,3];

u='x*y*z';

程序调用：s=example(node,elem,u)

程序输出：s =0.014433756729741

基金项目

贵州省科学技术基金项目《板振动问题的非协调有限元自适应算法》(项目合同编号：黔科合LH字[2014]7061号)。

附录

知网检索的两种方式：

1. 打开知网页面http://kns.cnki.net/kns/brief/result.aspx?dbPrefix=WWJD

下拉列表框选择：[ISSN]，输入期刊ISSN：2324-7991，即可查询

2. 打开知网首页http://cnki.net/

左侧“国际文献总库”进入，输入文章标题，即可查询

投稿请点击：http://www.hanspub.org/Submission.aspx

期刊邮箱：[e-mail]