二维定常对流扩散方程的重心插值配点法

Barycentric Interpolation Collocation Method for the Two-Dimensional Steady Convection-Diffusion Equation

摘要: 利用重心插值配点法求解二维定常对流扩散方程。首先介绍了两种重心插值配点法，并给出微分矩阵。其次，离散二维定常对流扩散方程以及初边值条件，利用置换法和附加法处理边界条件。采用第二类Chebyshev节点和等距节点进行数值计算，比较了两种边界条件施加方法下两种重心插值法的数值算法。数值算例表明了重心插值配点法的高精度性。

Abstract: The two-dimensional steady convection-diffusion equation is solved by barycentric interpolation method. Firstly, the two barycentric interpolation collocation methods are introduced, and the differential matrices are given. Secondly, the two-dimensional steady convection-diffusion equation and initial boundary conditions are dispersed, and the boundary conditions are treated by substitution and addition. Numerical calculations are carried out by using the second type of Chebyshev node and the equidistant node. Numerical examples demonstrate that this barycentric interpolation collocation method has high accuracy.

1. 引言

对流扩散方程在流体运动、传热和天体物理学、海洋学、气象学以及许多其他领域[1]-[3]都有广泛的应用，它的研究对科学和工程领域都具有重要的意义。

近年来，对于对流扩散方程数值求解已有大量研究，很多研究方法已被提出。数值计算方法主要包括：有限差分法[4]、有限元法[5]、有限体积法[6]等。Tian和Dai [7]提出了定常对流扩散方程的高阶紧致指数有限差分方法，求解了线性和非线性问题，Liao和Zhu [8]建立了四阶紧致差分格式求解了一维非定常对流扩散方程，消除对流项，结合中心差分格式求解，一方面提高了数值精度，另外可以将该方法可以直接应用于更高维的问题。

重心有理插值配点法和重心Lagrange插值法是无网格法的一种。Berrut [9]系统提出了重心Lagrange插值，并且在文中说明了重心插值法的收敛性和重心插值本身无关，而是取决于多项式插值函数，只要在区间[−1, 1]上的函数足够光滑，那么多项式插值就能以很快的收敛速度收敛到所求函数。同年，Higham在文献[10]中研究了重心Lagrange插值法是向前稳定的，对插值多项式的公式进行的误差分析，并且说明了在选择恰当的插值节点时重心插值公式的高精度性。2007年王兆清[11]等人将Lagrange插值公式改写为重心型插值公式，提出了重心型插值公式。Mittal和Rohila [12]将Burgers方程和Fisher方程转化为常微分方程组，应用重心有理插值法结合SSP-RK方法进行求解。王磊磊[13]利用重心插值配点法求解了二维广义变系数Poisson方程，推导出了计算公式，给出数值算例验证了两种重心插值配点法的优越性。文献[14]基于重心Lagrange插值法和重心有理插值法提出了两种新的求积方法，对Volterra积分方程求解。

本文采用两种重心插值配点法求解二维定常对流扩散方程。第2节介绍二维定常对流扩散问题。第3节介绍两种重心插值配点法。第4节给出本文所求解方程的计算过程。第5节给出数值算例验证本文所提方法的高精度性。

2. 二维定常对流扩散方程

本节考虑如下二维定常对流扩散方程：

$-a\frac{{\partial }^{2}u}{\partial x^2}-b\frac{{\partial }^{2}u}{\partial y^2}+c\frac{\partial u}{\partial x}+d\frac{\partial u}{\partial y}=f\left(x,y\right),\left(x,y\right)\in \Omega$ (1)

其中，$a$ 和$b$ 为扩散项系数，$c$ 和$d$ 为对流项系数，$f\left(x,y\right)$ 为源项，该问题的定义域为$\Omega =\left[a,b\right]\times \left[c,d\right]$ ，并且未知函数$u\left(x,y\right)$ 满足如下的边界条件：

$u\left(a,y\right)={\phi }_1\left(y\right),u\left(b,y\right)={\phi }_2\left(y\right)$ (2)

$u\left(x,c\right)={\varphi }_1\left(x\right),u\left(x,d\right)={\varphi }_2\left(x\right)$ (3)

3. 重心插值配点法及微分矩阵

3.1. 重心Lagrange插值法

给定$n+1$ 个不同的插值节点$x_i$ 以及对应的函数值为$u_i=u\left(x_i\right),i=0,1,2,\cdots ,n$ ，那么Lagrange插值多项式为

$u\left(x\right)={\displaystyle \sum }_{i=0}^n{l}_i\left(x\right)u_i$ (4)

其中Lagrange插值基函数：

$l_i\left(x\right)=\frac{{\displaystyle \prod }_{k=0,k\ne i}^n\left(x-x_k\right)}{{\displaystyle \prod }_{k=0,k\ne i}^n\left(x_i-x_k\right)}.$ (5)

令

$l\left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_n\right).$ (6)

则公式(4)可以改写为：

$u\left(x\right)={\displaystyle \sum _{i=0}^n\left(\frac{l\left(x\right)}{x-x_i}{\displaystyle \prod _{k=0,k\ne i}^n\frac{1}{x_i-x_k}}\right)}{u}_i$ (7)

定义重心插值权

${\omega }_i=\frac{1}{{\displaystyle \prod }_{k=0,k\ne i}^n\left(x_i-x_k\right)},$ (8)

则(7)式改写为

$u\left(x\right)={\displaystyle \sum }_{i=0}^n\left(\frac{l\left(x\right)}{x-x_i}{\omega }_i\right)u_i={\displaystyle \sum }_{i=0}^n{L}_i\left(x\right)u_i,i=0,1,2,\cdots ,n,$ (9)

公式(9)中的插值基函数为$L_i\left(x\right)=l\left(x\right)\frac{{\omega }_i}{x-x_i},$ 公式(9)是Lagrange插值的另外一种表达形式。

若让公式(9)中的$u_i=1$ 恒成立，则有：

$1={\displaystyle \sum }_{i=0}^n{L}_i\left(x\right)=l\left(x\right){\displaystyle \sum }_{i=0}^n\frac{{\omega }_i}{x-x_i},$ (10)

结合公式(9)和(10)，即有重心Lagrange插值公式：

$u\left(x\right)=\frac{{\displaystyle \sum }_{i=0}^n\frac{{\omega }_i}{x-x_i}{u}_i}{{\displaystyle \sum }_{i=0}^n\frac{{\omega }_i}{x-x_i}}={\displaystyle \sum }_{i=0}^n{L}_i\left(x\right)u_i.$ (11)

其中$L_i\left(x\right)=\frac{\frac{{\omega }_i}{x-x_i}}{{\displaystyle \sum }_{i=0}^n\frac{{\omega }_i}{x-x_i}}$ 为重心Lagrange插值基函数。

3.2. 重心有理插值法

在区间$\left[a,b\right]$ 上存在$n+1$ 个不同的计算节点$x_i,i=0,1,\cdots ,n,$ 且其对应的函数值为$u\left(x_i\right)$ ，简化为$u_i$ 。

我们选择合适的参数$d,0\le d\le n,d\in N,$ 对于节点$i=0,1,2,\cdots ,n-d,$ 选择$\left(x_i,u_i\right),\left(x_{i+1},u_{i+1}\right),\cdots \left(x_{i+d},u_{i+d}\right)$ 这$d+1$ 节点来计算对应的Lagrange插值多项式，即得

$u_i\left(x\right)={\displaystyle \sum }_{k=i}^{i+d}{\displaystyle \prod }_{j=i,j\ne k}^{i+d}\frac{x-x_j}{x_k-x_j}{u}_k$ (12)

构造$u_i\left(x\right)$ 的权函数${\lambda }_i\left(x\right),{\lambda }_i\left(x\right)=\frac{{\left(-1\right)}^i}{\left(x-x_i\right)\cdots \left(x-x_{i+d}\right)}$ ；根据公式(12)和权函数${\lambda }_i\left(x\right)$ 构造出重心有理插值函数

$u\left(x\right)=\frac{{\displaystyle \sum }_{i=0}^{n-d}{\lambda }_i\left(x\right)u_i\left(x\right)}{{\displaystyle \sum }_{i=0}^{n-d}{\lambda }_i\left(x\right)}=\frac{{\displaystyle \sum }_{i=0}^{n-d}{\left(-1\right)}^i{\displaystyle \sum }_{k=i}^{i+d}\frac{1}{x-x_k}{\displaystyle \prod }_{j=i,j\ne k}^{i+d}\frac{1}{x_k-x_j}{u}_k}{{{\displaystyle \sum }}_{i=0}^{n-d}{\lambda }_i\left(x\right)}=\frac{{\displaystyle \sum }_{k=0}^n\frac{{\omega }_k}{x-x_k}{u}_k}{{{\displaystyle \sum }}_{i=0}^{n-d}{\lambda }_i\left(x\right)}.$ (13)

其中${\omega }_k={\displaystyle \sum }_{i\in J_k}{\left(-1\right)}^i{\displaystyle \prod }_{j=i,j\ne k}^{i+d}\frac{1}{x_k-x_j}$ 为插值权，$J_k=\left\{i\in I:k-d\le i\le k\right\},I=\left\{0,1,\cdots ,n\right\},$ 根据$1={\displaystyle \sum }_{k=i}^{i+d}{\displaystyle \prod }_{j=i,j\ne k}^{i+d}\frac{x-x_j}{x_k-x_j},$ 可以得到：

${\displaystyle \sum }_{i=0}^{n-d}{\lambda }_i\left(x\right)={\displaystyle \sum }_{k=0}^n\frac{{\omega }_k}{x-x_k},$ (14)

则可以得到重心有理插值公式：

$u\left(x\right)=\frac{{\displaystyle \sum }_{k=0}^n\frac{{\omega }_k}{x-x_k}{u}_k}{{\displaystyle \sum }_{k=0}^n\frac{{\omega }_k}{x-x_k}}={\displaystyle \sum }_{k=0}^n{L}_k\left(x\right)u_k.$ (15)

其中$L_k\left(x\right)=\frac{\frac{{\omega }_k}{x-x_k}}{{\displaystyle \sum }_{k=0}^n\frac{{\omega }_k}{x-x_k}}$ 为重心有理插值基函数。

3.3. 微分矩阵

两种重心插值配点法在形式上是一致的，根本不同在于权函数的不同，使用重心有理插值方法进行计算时，需要选取合适的$d$ 。

对于区间$\left[a,b\right]$ 上的节点$$ $a=x_1<x_2<\cdots <x_n=b,$ 函数$u\left(x,y\right)$ 在节点处的函数值为$u_j=u\left(x_j\right),j=1,2,\cdots ,n,$ 则函数$u\left(x\right)$ 的重心型插值函数为

$u\left(x\right)={\displaystyle \sum }_{j=1}^n{L}_j\left(x\right)u_j$ (16)

对公式(16)在节点$\left\{x_i,i=1,2,\cdots ,n\right\}$ 处求$k$ 阶导数，所得结果为

$u^{\left(k\right)}\left(x_i\right)={\displaystyle \sum }_{j=1}^n{L}_j^{\left(k\right)}\left(x_i\right)u_j,i=1,2,\cdots ,n$

将其写成简洁的微分矩阵如下

$u^{\left(k\right)}=C^{\left(k\right)}u,k=1,2,\cdots$ (17)

其中$u^{\left(k\right)}={\left[u^{\left(k\right)}\left(x_1\right),u^{\left(k\right)}\left(x_2\right),\cdots u^{\left(k\right)}\left(x_n\right)\right]}^T,u={\left[u_1,u_2,\cdots ,u_n\right]}^T,$

$C^{\left(k\right)}=\left(\begin{array}{cccc}{L}_1^{\left(k\right)}\left(x_1\right)& L_2^{\left(k\right)}\left(x_1\right)& \cdots & L_n^{\left(k\right)}\left(x_1\right)\\ L_1^{\left(k\right)}\left(x_2\right)& L_2^{\left(k\right)}\left(x_2\right)& \cdots & L_n^{\left(k\right)}\left(x_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ L_1^{\left(k\right)}\left(x_n\right)& L_2^{\left(k\right)}\left(x_n\right)& \cdots & L_n^{\left(k\right)}\left(x_n\right)\end{array}\right).$ (18)

当$k=0$ 时，$C^{\left(0\right)}=I_n,$ $I_n$ 是$n$ 阶单位矩阵。

当$k=1$ 时插值基函数的$1$ 阶导数写成统一形式如下：

$L_j^{\text{'}}\left(x_i\right)=\{\begin{array}{r}\frac{\frac{{\omega }_j}{{\omega }_i}}{x_i-x_j},i\ne j\\ -{\displaystyle \sum }_{j=1,j\ne i}^n{L}_j^{\text{'}}\left(x_i\right),i=j\end{array}$ (19)

当$k\ge 2$ 时插值基函数的$k$ 阶导数写成统一形式如下：

$L_j^{\left(k\right)}\left(x_i\right)=\{\begin{array}{r}k\left(L_i^{\left(k-1\right)}\left(x_i\right)L_j^{\text{'}}\left(x_i\right)-\frac{L_j^{k-1}\left(x_i\right)}{x_i-x_j}\right),i\ne j\\ -{\displaystyle \sum }_{j=1,j\ne i}^n{L}_j^{\left(k\right)}\left(x_i\right),i=j\end{array}k=2,3,\cdots ,$ (20)

以上关于微分矩阵的详细推导过程见文献[15]中第1.6节。

3.4. 插值节点

由于两种插值公式中的插值权均和插值节点的选取有关，我们选取特殊的两种插值节点，即等距插值节点和Chebyshev插值节点进行计算，由于等距节点在Matlab中可以直接利用linspace指令生成，因此下面只给出Chebyshev节点的介绍。

(1) 第一类Chebyshev节点

$x_j=\cos \frac{\left(2j+1\right)\pi }{2n+2},j=0,1,2,\cdots ,n$

(2) 第二类Chebyshev节点

$x_j=\cos \frac{i\pi }{n},j=0,1,2,\cdots ,n$

需要注意的是，Chebyshev节点是定义在区间$\left[-1,1\right]$ 内的节点族，对于一般计算区间$\left[a,b\right]$ ，可以利用区间坐标变换公式$\overline{x}=\frac{\left(b-a\right)x}{2}+\frac{b+a}{2}$ 转换计算即可。

4. 重心插值配点法求解二维定常对流扩散方程

4.1. 双变量重心插值法的介绍

在此之前，我们先介绍双变量的重心型插值配点法。

关于变量$x$ 、$y$ 的未知函数$u\left(x,y\right)$ ，我们将定义区间$\Omega =\left[a,b\right]\times \left[c,d\right]$ 在$x$ 方向和$y$ 方向分别离散为$m$ 、$n$ 个计算节点，即$a=x_1<x_2<\cdots <x_m=b$ 和$c=y_1<y_2<\cdots <y_n=d$ ，那么在整个计算区域上便有$m\times n$ 个张量积型插值节点$\left\{\left(x_i,y_j\right),i=1,2,\cdots ,m;j=1,2,\cdots ,n\right\}$

那么未知函数$u\left(x,y\right)$ 的重心Lagrange插值公式可以表示成

$u\left(x,y\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x\right)M_j\left(y\right)u_{ij}}},$ (21)

其中$u_{ij}=u\left(x_i,y_j\right),L_i\left(x\right)、M_j\left(y\right)$ 分别代表方向$x$ 和$y$ 方向上的插值基函数，并且具体形式为：

$L_i\left(x\right)=\frac{\frac{{\omega }_i}{x-x_i}}{{\displaystyle \sum _{i=1}^m\frac{{\omega }_i}{x-x_i}}},M_j\left(\text{y}\right)=\frac{\frac{{\nu }_j}{y-y_j}}{{\displaystyle \sum _{j=1}^n\frac{{\nu }_j}{y-y_j}}},i=1,2,\cdots ,m;j=1,2,\cdots ,n$ (22)

$x$ 和$y$ 方向上的插值权重分别为：

${\omega }_i=\frac{1}{{\displaystyle \prod _{k=1,k\ne i}^m\left(x_i-x_k\right)}},{\nu }_j=\frac{1}{{\displaystyle \prod _{k=1,k\ne j}^n\left(y_j-y_k\right)}}$ (23)

为简化离散矩阵，接下来，对公式(21)求$l+r$ 阶偏导数有：

$\frac{{\partial }^{l+r}u\left(x,y\right)}{\partial x^l\partial y^r}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i^{\left(l\right)}\left(x\right)M_j^{\left(r\right)}\left(y\right)u_{ij}}}$ (24)

那么(24)式在插值节点$\left(x_j,y_k\right),j=1,2,\cdots ,m,k=1,2,\cdots ,n$ 处的$l+r$ 阶偏导数为：

$\frac{{\partial }^{l+r}u\left(x_j,y_k\right)}{\partial x^l\partial y^r}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i^{\left(l\right)}\left(x_j\right)M_j^{\left(r\right)}\left(y_k\right)u_{ij}}}$ (25)

写成矩阵的形式：

$u^{\left(l+r\right)}=C_x^{\left(l\right)}{D}_y^{\left(r\right)}u$ (26)

其中，$C_x^{\left(l\right)}$ 是$x$ 方向上的$l$ 阶微分矩阵，$D_y^{\left(r\right)}$ 是$y$ 方向上的$r$ 阶微分矩阵。

4.2. 重心Lagrange插值法的离散过程

将式(21)代入方程(1)中，得到

$\left[\begin{array}{l}-a{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L^{″}}_i\left(x\right)M_j\left(y\right)}}-b{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i(x){M^{″}}_j\left(y\right)}}+\\ c{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L^{\prime }}_i\left(x\right)M_j\left(y\right)}}+d{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x\right){M^{\prime }}_j\left(y\right)}}\end{array}\right]u_{ij}=f\left(x,y\right)$ (27)

令(27)在节点$\left(x_j,y_k\right),j=1,2,\cdots ,m;k=1,2,\cdots ,n$ 处均成立，将这些节点代入上式有：

$\left[\begin{array}{l}-a{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L^{″}}_i\left(x_j\right)M_j\left(y_k\right)}}-b{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x_j\right){M^{″}}_j\left(y_k\right)}}+\\ c{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L^{\prime }}_i\left(x_j\right)M_j\left(y_k\right)}}+d{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x_j\right){M^{\prime }}_j\left(y_k\right)}}\end{array}\right]u_{ij}=f\left(x_j,y_k\right)$ (28)

由插值基函数的性质以及微分矩阵的相关知识：

$L_i(x_j)={\delta }_{ij},{L^{\prime }}_i\left(x_j\right)=D_{ij}^{\left(1\right)},{L^{″}}_i\left(x_j\right)=D_{ij}^{\left(2\right)},M_j\left(y_k\right)={\delta }_{jk},{M^{\prime }}_j\left(y_k\right)=C_{ij}^{\left(1\right)},{M^{″}}_j\left(y_k\right)=C_{ij}^{\left(2\right)}.$

那么(28)式变成

$\left[-a\left(D^{\left(2\right)}\otimes I_n\right)-b\left(I_m\otimes C^{\left(2\right)}\right)+c\left(D^{\left(1\right)}\otimes I_n\right)+d\left(I_m\otimes C^{\left(1\right)}\right)\right]\left(\begin{array}{c}{u}_{1j}\\ ⋮\\ u_{mj}\end{array}\right)=\left(\begin{array}{c}{f}_{1j}\\ ⋮\\ f_{mj}\end{array}\right)$ (29)

令

$U={\left[u_{11},\cdots ,u_{1n},u_{21},\cdots ,u_{2n},\cdots ,u_{m1},\cdots ,u_{mn}\right]}^T$

$F={\left[f_{11},\cdots ,f_{1n},f_{21},\cdots ,f_{2n},\cdots ,f_{m1},\cdots ,f_{mn}\right]}^T$

那么(29)式变成

$\left[-a\left(D^{\left(2\right)}\otimes I_n\right)-b\left(I_m\otimes C^{\left(2\right)}\right)+c\left(D^{\left(1\right)}\otimes I_n\right)+d\left(I_m\otimes C^{\left(1\right)}\right)\right]U=F$ (30)

即

$LU=F$ (31)

其中

$L=-a\left(D^{\left(2\right)}\otimes I_n\right)-b\left(I_m\otimes C^{\left(2\right)}\right)+c\left(D^{\left(1\right)}\otimes I_n\right)+d\left(I_m\otimes C^{\left(1\right)}\right)$

接下来对边界条件进行离散：

利用重心Lagrange插值法对(2)进行近似，可以得到离散表达式

$u\left(a,y\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(a\right)M_j\left(y\right)}}{u}_{ij}={\phi }_1\left(y\right)$ (32)

$u\left(b,y\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(b\right)M_j\left(y\right)}}{u}_{ij}={\phi }_2\left(y\right)$ (33)

同理(3)式的离散表达式为：

$u\left(x,c\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x\right)M_j\left(c\right)}}{u}_{ij}={\varphi }_1\left(x\right)$ (34)

$u\left(x,d\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{L}_i\left(x\right)M_j\left(d\right)}}{u}_{ij}={\varphi }_2\left(x\right)$ (35)

将(32)-(35)写成矩阵的形式为：

$\left(I_{m1}\otimes I_n\right)U={\phi }_1,\left(I_{mm}\otimes I_n\right)U={\phi }_2$ (36)

$\left(I_m\otimes I_{n1}\right)U={\varphi }_1,\left(I_m\otimes I_{nn}\right)U={\varphi }_2$ (37)

其中，$I_m$ 和$I_n$ 分别表示$m$ 阶单位矩阵和$n$ 阶单位矩阵，$I_{m1}$ 和$I_{mm}$ 分别表示$m$ 阶单位矩阵的第1行和第$m$ 行，$I_{n1}$ 和$I_{nn}$ 分别表示$n$ 阶单位矩阵的第1行和第$n$ 行。

附加法施加边界条件：没有施加边界条件时(31)式是$m\times n$ 维矩阵，将(36)式中左边的$\left(I_{m1}\otimes I_n\right)$ 各行直接附加到(31)式等号的左边，等号右边直接代入，形成$m\times n+2m$ 维矩阵方程。同理，对其它式子也进行同样操作。

置换法施加边界条件：对于$x$ 方向上的边界条件，将(37)式中左边的$\left(I_m\otimes I_{n1}\right)$ 各行式子对应置换(31)式等号的左边矩阵中的第$\left\{1,\cdots ,n\right\}$ 行，将式子中$\left(I_m\otimes I_{nn}\right)$ 中各行对应置换(31)式等号的左边矩阵的$\left\{\left(m-1\right)\times n+1:m\times n\right\}$ 行，对应地，等式右边也进行置换。对于$y$ 方向上边界条件的置换也类似。

至此，我们完成了利用重心Lagrange插值法离散二维定常对流扩散方程的过程以及边界条件的处理。

4.3. 重心有理插值法的离散过程

同理，将区间定义为$\Omega =\left[a,b\right]\times \left[c,d\right]$ ，结合式(21)，我们可以得到关于变量$x,y$ 的未知函数$\overline{u}\left(x,y\right)$ 的重心有理插值公式为

$\overline{u}\left(x,y\right)={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{\overline{L}}_i\left(x\right){\overline{M}}_j\left(y\right)u_{ij}}}$ (38)

其中$u_{ij}=\overline{u}\left(x_i,y_j\right),{\overline{L}}_i\left(x\right)、{\overline{M}}_j\left(y\right)$ 分别代表$x$ 方向和$y$ 方向上的插值基函数，并且具体形式和重心Lagrange插值法的类似，${\overline{\omega }}_i$ 和${\overline{v}}_j$ 分别是$x$ 和$y$ 方向上的插值权重。

将(38)代入方程(1)中，经过化简得到如下计算格式：

$\left[-a\left({\overline{D}}^{\left(2\right)}\otimes I_n\right)-b\left(I_m\otimes {\overline{C}}^{\left(2\right)}\right)+c\left({\overline{D}}^{\left(1\right)}\otimes I_n\right)+d\left(I_m\otimes {\overline{C}}^{\left(1\right)}\right)\right]U=F$ (39)

即

$\overline{L}U=F$ (40)

其中

$\overline{L}=-a\left({\overline{D}}^{\left(2\right)}\otimes I_n\right)-b\left(I_m\otimes {\overline{C}}^{\left(2\right)}\right)+c\left({\overline{D}}^{\left(1\right)}\otimes I_n\right)+d\left(I_m\otimes {\overline{C}}^{\left(1\right)}\right)$

$U={\left[u_{11},\cdots ,u_{1n},u_{21},\cdots ,u_{2n},\cdots ,u_{m1},\cdots ,u_{mn}\right]}^T$

$F={\left[f_{11},\cdots ,f_{1n},f_{21},\cdots ,f_{2n},\cdots ,f_{m1},\cdots ,f_{mn}\right]}^T$

同重心Lagrange插值法，对边界条件使用重心有理插值法进行离散：

可以得到二维定常对流扩散方程$x$ 的边界条件以及y的边界条件的最终离散矩阵形式为：

$\left({\overline{I}}_{m1}\otimes I_n\right)U={\phi }_1,\left({\overline{I}}_{mm}\otimes I_n\right)U={\phi }_2$ (41)

$\left(I_m\otimes {\overline{I}}_{n1}\right)U={\varphi }_1,\left(I_m\otimes {\overline{I}}_{nn}\right)U={\varphi }_2$ (42)

其中，${\overline{I}}_m$ 和${\overline{I}}_n$ 分别表示$m$ 阶单位矩阵和$n$ 阶单位矩阵，${\overline{I}}_{m1}$ 和${\overline{I}}_{mm}$ 分别表示$m$ 阶单位矩阵的第1行和第$m$ 行，${\overline{I}}_{n1}$ 和${\overline{I}}_{nn}$ 分别表示$n$ 阶单位矩阵的第1行和第$n$ 行。

附加法和置换法施加边界条件同重心Lagrange插值法离散部分的介绍，这里不再赘述。

5. 数值算例

在进行数值计算之前，我们首先给出绝对误差和相对误差的定义：

(1) 绝对误差：

$Er=\Vert u_c-u_e\Vert$

(2) 相对误差：

$\mathrm{Re}\text{r}=\frac{Er}{\Vert u_e\Vert }$

其中$\Vert •\Vert$ 表示向量的范数，$u_e$ 表示精确解，$u_c$ 表示数值解。

下面来验证重心型插值配点法求解二维定常对流扩散方程的高精度性，我们给出如下算例。

算例 考虑二维定常对流扩散方程

$-2\frac{\partial u}{\partial x}-2\frac{\partial u}{\partial y}-\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=f\left(x,y\right)$

其中$\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right]$ ，

$f\left(x,y\right)=e^{\left(x+y\right)}\left[\left(6-5{\pi }^2\right)\sin \left(\pi x\right)\sin \left(2\pi y\right)+4\pi \cos \left(\pi x\right)\sin \left(2\pi y\right)+8\pi \sin \left(\pi x\right)\cos \left(2\pi y\right)\right]$

该方程具有解析解

$u\left(x,y\right)=e^{\left(x+y\right)}\sin \left(\pi x\right)\sin \left(2\pi y\right)$ 。

边界条件由解析解确定。

取插值节点数为$19\times 19$ ，$d=17$ 。表1给出两种重心插值法在两种不同边界条件施加方法下的计算误差，结果显示，使用置换法施加边界条件时，两种重心插值法的计算效果均比较好；使用附加法施加边界条件时，重心有理插值法和重心有理插值法的计算误差都比置换法下的计算误差大，因此在本算例中我们选取置换法进行研究。

Table 1. Errors in different boundary condition imposition methods for two centre of gravity interpolation methods

表1. 两种重心插值法不同边界条件施加方法的误差

施加方法	重心Lagrange插值		重心有理插值
	Er	Rer	Er	Rer
置换法	4.7323e-12	2.3450e-13	6.9204e-12	3.4294e-13
附加法	1.9005e-09	9.4179e-11	1.8991e-09	9.4110e-11

将上述结果与文献[16]中的数值计算方法积分方程法和有限体积法分别进行比较，对比结果见表2。我们发现两种重心插值法求解对流扩散方程时在较小节点数下就能够取得较小的计算误差。而文献[16]中所提方法在节点数达到$100\times 100$ 时，两种计算数值计算方法的计算误差数量级分别为${10}^{-5}$ 和${10}^{-4}$ ，由此可以表明重心插值法的高精度性。但是重心插值法在处理多个边界条件或者复杂边界条件在编程上的处理比较复杂。

Table 2. Comparison of computational errors of several numerical calculation methods

表2. 几种数值计算方法的计算误差的比较

计算方法	Er	Rer
积分方程法	1.7002e-3	3.2920e-4
有限体积法	1.8961e-2	4.5692e-3
重心Lagrange插值法	4.7323e-12	2.3450e-13
重心有理插值法	6.9204e-12	3.4294e-13

表3给出两种重心插值法在不同插值节点类型下随节点数的增加计算误差的变化，计算结果显示，在第二类Chebyshev节点下，随着节点数的增加，重心Lagrange插值法的计算误差变化幅度较小，误差数量级保持在${10}^{-11}$ ，整体上保持不错的计算精度，重心有理插值法的计算误差随节点数的增加逐渐变大，在节点数达到$35\times 35$ 时，误差数量级保持在${10}^{-10}$ ；在等距节点下，重心Lagrange插值法的计算误差随着节点数的增加越来越大，当节点数达到$25\times 25$ 时，数量级变成${10}^{-1}$ ，节点数再增加时，重心Lagrange插值法所得误差变得非常大，数值算法不稳定；而重心有理插值法的计算误差随节点数的增加先变小再变大，误差数量级最后保持在${10}^{-5}$ ，结果比重心Lagrange插值法下的结果要好，因此，两种重心插值法在第二类Chebyshev节点下均有良好的计算精度，而在等距节点下，重心Lagrange插值法更适用于较少数量的插值节点，重心有理插值法的适用性更好。

Table 3. Computational errors of two interpolation methods under two types of interpolation nodes

表3. 两类插值节点下两种插值法的计算误差

节点数	等距节点		第二类Chebyshev节点
	重心Lagrange	重心有理	重心Lagrange	重心有理
$21\times 21$	6.6145e-05	1.1592e-05	9.4067e-12	1.3138e-11
$23\times 23$	4.8583e-03	9.2491e-06	1.0495e-11	1.0644e-11
$25\times 25$	2.2330e-01	1.5676e-05	1.4978e-11	2.7430e-11
$27\times 27$	5.6331e-01	2.5067e-05	2.1299e-11	2.3914e-11
$29\times 29$	1.4909e+02	2.0283e-05	2.5357e-11	4.5984e-11
$31\times 31$	4.0335e+03	2.1826e-05	3.4243e-11	8.7449e-11
$33\times 33$	9.6827e+03	3.2130e-05	3.7722e-11	1.1507e-10
$35\times 35$	5.3974e+04	7.3608e-05	7.8012e-11	2.3820e-10

根据前面的分析，给出两种重心插值法在第二类Chebyshev节点下的数值解如图1，结果显示两种插值法都具有良好的数值拟合性。

(a) 重心Lagrange插值法的数值解图 (b) 重心有理插值法的数值解图

Figure 1. Numerical solution plots of the two interpolation methods under Chebyshev node

图1. Chebyshev节点下两种插值方法的数值解图

下面给出第二类Chebyshev节点下两种重心插值法对应的误差分布图，见图2，重心Lagrange插值法在区域边界处的震荡要比重心有理插值法下的震荡大一些，两种重心插值法在中心区域处均比较稳定。

(a) 重心Lagrange插值法的误差图 (b) 重心有理插值法的误差图

Figure 2. Error distribution of two interpolation methods under Chebyshev node

图2. Chebyshev节点下两种插值法的误差分布

6. 结论

本文介绍求解二维定常对流扩散方程的重心Lagrange插值配点法和重心有理插值配点法，得到定常对流扩散方程离散化的矩阵格式。数值算例验证了两种重心插值方法在Chebyshev节点下的数值解和精确解均具有较好的吻合度。在等距节点下，重心有理插值法的计算精度要比重心Lagrange插值法的计算精度高。

基金项目

陕西省科技厅重点研发一般资助项目(2024M722604)；西安市科技局高校院所科技人员服务企业项目(24GXFW0038)。

NOTES

^*通讯作者。