重心插值配点法求解二维波动方程

Barycentric Interpolation Collocation Method for the Two-Dimensional Wave Equation

摘要

利用重心有理插值配点法求解波动方程。首先，介绍重心有理插值配点法并且给出微分矩阵；其次，采用该方法离散二维波动方程和初边值条件，利用置换法处理边界条件，获得方程的最终矩阵形式。选取第二类Chebyshev节点和等距节点进行数值计算，将重心有理插值配点法得到的数值精度与重心Lagrange插值配点法得到的数值精度进行了比较，数值算例表明，两种方法在Chebyshev节点下均能保持高精度性以及稳定性，而在等距节点下重心有理插值法的数值计算效果更好。

Abstract

The barycentric rational interpolation collocation method is used to solve the wave equation. Firstly, the method is introduced and the differential matrix is given. Secondly, the method is used to discretize the two-dimensional wave equation and the initial boundary conditions, and the replacement method is used to deal with the boundary conditions to obtain the final matrix form of the equation. The second kind of Chebyshev nodes and the equidistant nodes are chosen for numerical computation. The numerical accuracy of the method in this study is compared with that by the barycentric Lagrange interpolation, and the numerical example demonstrate that both methods maintain high accuracy as well as stability under Chebyshev nodes, while the barycentric rational interpolation method works better numerically under equidistant nodes.

Graphical abstract

关键词

波动方程 / 重心有理插值 / Chebyshev节点 / 等距节点

Key words

wave equation / barycentric rational interpolation / Chebyshev node / equidistant node

波动方程是一种重要的双曲型微分方程，常用来描述自然界内各种波现象，例如声波、光波、水波^[1-4]，其被广泛应用于弹性力学、地球物理、电磁学以及气象学等诸多自然科学领域^[5-6]。关于波动方程数值求解的研究一直以来广受学者的关注。

波动方程的主要数值计算方法有：有限差分法^[7-8]、有限元法^[9]、有限体积法^[10]和谱方法^[11]。文献[12]开发了一种有限差分法结合交替方向隐式方案，求解了具有分数阻尼的二维波动方程。文献[13]采用有限元方法进行空间离散，时间导数采用有限差分法近似，解决了扩散黏性波动方程的一般初始边值问题。

无网格法是一种高精度算法，常用来处理微分方程的初边值问题。其中重心有理插值配点法就是无网格法的一种。王兆清等^[14]将Lagrange插值公式改写为重心型插值公式，提出一维重心型插值公式。马燕等^[15]提出了(1+1)维波动方程的重心插值配点法，采用重心有理插值函数离散波动方程。宋灵宇等^[16]采用重心插值配点法求解了二维Sobolev方程，取不同插值节点进行计算，实验结果表明重心有理插值法在不同插值节点下都能够获得稳定性较好的数值解。屈金铮等^[17]利用重心Lagrange插值配点法求解非线性方程，采用三种迭代格式求解了伪抛物方程。袁洪旺等^[18]提出了求解高维波动方程的重心Lagrange插值配点法，分别使用置换法和附加法施加初边值条件求解波动方程。文献[19]介绍了重心有理插值法求解(1+1)维Burgers方程，并且证明了Burgers方程的重心有理插值配点法的收敛速度。文献[20]选择重心有理插值配点法求解SG方程，构造直接线性化等迭代格式，比较了不同插值节点下各迭代格式的结果，给出数值算例验证了方法的有效性。

本文采用重心有理插值配点法求解二维波动方程。第1节介绍本文所要求解的二维波动方程。第2节介绍重心有理插值配点法。第3节给出波动方程的计算过程，包括波动方程的离散格式以及初边值条件的离散。第4节给出数值算例验证本文所提方法的高精度和高效性。第5节给出了结论。

1 二维波动方程

本节主要考虑如下形式的二维波动问题：

\left\{\begin{array}{l}{u}_{tt}=c^{\mathrm{2}}\left(u_{xx}+u_{yy}\right)+f\left(x,y,t\right)\\ u\left(x,y,\mathrm{0}\right)=\varphi \left(x,y\right),\left(x,y\right)\in \Omega \\ u_t\left(x,y,\mathrm{0}\right)=\psi \left(x,y\right),\left(x,y\right)\in \Omega \\ u\left(x,y,t\right)=g\left(x,y,t\right),\left(x,y\right)\in \Gamma \end{array}\right.(1)

其中：(x,y,t)\in \Omega \times (\mathrm{0},T]，\Omega =[a,b]\times [c,d]，\Gamma为\Omega的边界；c^{\mathrm{2}}是波动系数且为已知常数，u(x,y,t)为未知函数，f(x,y,t)为非齐次项；\varphi (x,y,t)，\psi (x,y,t)和g\left(x,y,t\right)均为已知函数。

2 重心有理插值配点法及微分矩阵

2.1 重心有理插值

设在区间[a,b]上有n+\mathrm{1}个不同的计算节点x_i，i=\mathrm{0},\mathrm{1},\cdots ,n且其对应的函数值为u(x_i)，记u_i。

选择适当的参数d,\mathrm{0}\le d\le n,d\in N，对d+\mathrm{1}个节点(x_i,u_i),(x_{i+\mathrm{1}},u_{i+\mathrm{1}}),\cdots ,(x_{i+d},u_{i+d})计算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(2)

构造多项式u_i(x)的权函数

{\lambda }_i(x)=\frac{{\left(-\mathrm{1}\right)}^i}{\left(x-x_i\right)\cdots \left(x-x_{i+d}\right)}

根据式(2)和权函数{\lambda }_i(x)构造重心有理插值函数：

\begin{array}{l}u\left(x\right)=\frac{{\displaystyle \sum _{i=\mathrm{0}}^{n-d}}{\lambda }_i\left(x\right)u_i\left(x\right)}{{\displaystyle \sum _{i=\mathrm{0}}^{n-d}}{\lambda }_i\left(x\right)}\\ =\frac{{\displaystyle \sum _{i=\mathrm{0}}^{n-d}}{\left(-\mathrm{1}\right)}^i{\displaystyle \sum _{k=i}^{i+d}}\frac{\mathrm{1}}{x-x_k}{\displaystyle \prod _{j=i,j\ne k}^{i+d}}\frac{\mathrm{1}}{x_k-x_j}{u}_k}{{\displaystyle \sum _{i=\mathrm{0}}^{n-d}}{\lambda }_i\left(x\right)}\\ =\frac{{\displaystyle \sum _{k=\mathrm{0}}^n}\frac{{\omega }_k}{x-x_k}{u}_k}{{\displaystyle \sum _{i=\mathrm{0}}^{n-d}}{\lambda }_i\left(x\right)}\end{array}(3)

其中：{\omega }_k=\sum _{i\in J_k}{\left(-\mathrm{1}\right)}^i{\displaystyle \prod _{j=i,j\ne k}^{i+d}}\frac{\mathrm{1}}{x_k-x_j}为重心有理插值权,J_k=\left\{i\in I:k-d\le i\le k\right\},I=\left\{\mathrm{0},\mathrm{1},\cdots ,n\right\}，根据\mathrm{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=\mathrm{0}}^{n-d}}{\lambda }_i\left(x\right)={\displaystyle \sum _{k=\mathrm{0}}^n}\frac{{\omega }_k}{x-x_k}(4)

因此，由式(3)、式(4)可以得到重心有理插值公式

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

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

2.2 微分矩阵

使用重心有理插值方法进行计算时，d会影响计算的精度，因此需要选取合适的d值。

对式(5)在节点\left\{x_i,i=\mathrm{1},\mathrm{2},\cdots ,n\right\}处求导：

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

将其写成简洁的微分矩阵

{\mathit{u}}^{\left(p\right)}={\mathit{C}}^{\left(p\right)}\mathit{u},p=\mathrm{1},\mathrm{2},\cdots(6)

其中：

\begin{array}{l}{\mathit{u}}^{\left(p\right)}={\left[{\mathit{u}}^{\left(p\right)}\left(x_{\mathrm{1}}\right),{\mathit{u}}^{\left(p\right)}\left(x_{\mathrm{2}}\right),\cdots ,{\mathit{u}}^{\left(p\right)}\left(x_n\right)\right]}^{\mathrm{T}}\\ \mathit{u}={\left[u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_n\right]}^{\mathrm{T}}\\ {\mathit{C}}^{\left(p\right)}=\left[\begin{array}{cccc}{L}_{\mathrm{1}}^{\left(p\right)}\left(x_{\mathrm{1}}\right)& L_{\mathrm{2}}^{\left(p\right)}\left(x_{\mathrm{1}}\right)& \cdots & L_n^{\left(p\right)}\left(x_{\mathrm{1}}\right)\\ L_{\mathrm{1}}^{\left(p\right)}\left(x_{\mathrm{2}}\right)& L_{\mathrm{2}}^{\left(p\right)}\left(x_{\mathrm{2}}\right)& \cdots & L_n^{\left(p\right)}\left(x_{\mathrm{2}}\right)\\ ⋮& ⋮& & ⋮\\ L_{\mathrm{1}}^{\left(p\right)}\left(x_n\right)& L_{\mathrm{2}}^{\left(p\right)}\left(x_n\right)& \cdots & L_n^{\left(p\right)}\left(x_n\right)\end{array}\right]\end{array}(7)

当p=\mathrm{0}时，{\mathit{C}}^{\left(\mathrm{0}\right)}={\mathit{I}}_n,{\mathit{I}}_n是n阶单位矩阵。

当p=\mathrm{1}时，插值基函数的1阶导数写成统一形式如下：

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

当p\ge \mathrm{2}时，插值基函数的p阶导数写成统一形式如下：

L_j^{\left(p\right)}\left(x_i\right)=\left\{\begin{array}{l}p\left[L_i^{\left(p-\mathrm{1}\right)}\left(x_i\right)L_j^{\mathrm{\text{'}}}\left(x_i\right)-\frac{L_j^{p-\mathrm{1}}\left(x_i\right)}{x_i-x_j}\right],i\ne j\\ -{\displaystyle \sum _{j=\mathrm{1},j\ne i}^n}{L}_j^{\left(p\right)}\left(x_i\right),i=j\end{array}\right.(9)

3 波动方程的计算

将重心有理插值配点法运用到二维波动方程的求解中。

考虑函数u\left(x,y,t\right),\left(x,y,t\right)\in \Omega \times \left[\mathrm{0},T\right]，其中\Omega =\left[a,b\right]\times \left[c,d\right]，将\Omega离散为m\times n个张量型插

值节点\left\{{\Omega }_{ij}=\left(x_i,y_j\right),i=\mathrm{1},\mathrm{2},\cdots ,m;j=\mathrm{1},\mathrm{2},\cdots ,n\right\}，将

\left[\mathrm{0},T\right]离散为l个插值节点\left\{t_k,k=\mathrm{1},\mathrm{2},\cdots ,l\right\}，最终在\Omega \times \left[\mathrm{0},T\right]内得到m\times n\times l个张量型插值节点：

\left\{\begin{array}{l}\left(x_i,y_j,t_k\right),i=\mathrm{1},\mathrm{2},\cdots ,m;\\ j=\mathrm{1},\mathrm{2},\cdots ,n;t=\mathrm{1},\mathrm{2},\cdots ,l\end{array}\right\}

将函数u\left(x,y,t\right)在插值节点\left(x_i,y_j,t_k\right)处的函数值定义为u\left(x_i,y_j,t_k\right)=u_{ijk}，则u\left(x,y,t\right)在节点\left(x_i,y_j,t_k\right)上的重心有理插值表达式为

u\left(x,y,t\right)={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^l}{L}_i\left(x\right)M_j\left(y\right)T_k\left(t\right)u_{ijk}(10)

其中

\begin{array}{l}{L}_i\left(x\right)=\frac{\frac{{\omega }_i}{x-x_i}}{{\displaystyle \sum _{i=\mathrm{1}}^m}\frac{{\omega }_i}{x-x_i}}\\ M_j\left(y\right)=\frac{\frac{{\nu }_j}{y-y_j}}{{\displaystyle \sum _{j=\mathrm{1}}^n}\frac{{\nu }_j}{y-y_j}}\\ T_k\left(t\right)=\frac{\frac{{\lambda }_k}{t-t_k}}{{\displaystyle \sum _{k=\mathrm{1}}^l}\frac{{\lambda }_k}{t-t_k}}\end{array}(11)

分别是x，y，t方向上的重心有理插值基函数。{\omega }_i,

{\nu }_j,{\lambda }_k分别为对应的重心有理插值权函数：

\begin{array}{l}{\omega }_i=\frac{\mathrm{1}}{{\displaystyle \prod _{k=\mathrm{1},k\ne i}^m}\left(x_i-x_k\right)}\\ {\nu }_j=\frac{\mathrm{1}}{{\displaystyle \prod _{k=\mathrm{1},k\ne j}^n}\left(y_j-y_k\right)}\\ {\lambda }_k=\frac{\mathrm{1}}{{\displaystyle \prod _{i=\mathrm{1},i\ne k}^l}\left(t_k-t_j\right)}\end{array}(12)

进一步，式(10)的\left(g_{\mathrm{1}}+g_{\mathrm{2}}+g_{\mathrm{3}}\right)阶偏导数可以表示为：

\begin{array}{c}{u}^{\left(g_{\mathrm{1}}+g_{\mathrm{2}}+g_{\mathrm{3}}\right)}\left(x,y,t\right)={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^l}{L_i}^{\left(g_{\mathrm{1}}\right)}\left(x\right)\\ {M_j}^{\left(g_{\mathrm{2}}\right)}\left(y\right){T_k}^{\left(g_{\mathrm{3}}\right)}\left(t\right)u_{ijk}\end{array}(13)

则式(13)在插值节点\left\{\begin{array}{l}\left(x_r,y_s,t_p\right),r=\mathrm{1},\mathrm{2},\cdots ,m;\\ p=\mathrm{1},\mathrm{2},\cdots ,l\end{array}\right.\left.\begin{array}{l}s=\mathrm{1},\mathrm{2},\cdots ,n\\ \\ \end{array}\right\}处的函数值可以离散为

\begin{array}{l}{u}^{\left(g_{\mathrm{1}}+g_{\mathrm{2}}+g_{\mathrm{3}}\right)}\left(x_r,y_s,t_p\right)\\ ={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^l}{L_i}^{\left(g_{\mathrm{1}}\right)}\left(x_r\right){M_j}^{\left(g_{\mathrm{2}}\right)}\left(y_s\right){T_k}^{\left(g_{\mathrm{3}}\right)}\left(t_p\right)u_{ijk}\end{array}(14)

写成微分矩阵的形式为

{\mathit{u}}^{\left(g_{\mathrm{1}}+g_{\mathrm{2}}+g_{\mathrm{3}}\right)}=\left[{\mathit{C}}^{\left(g_{\mathrm{1}}\mathrm{00}\right)}\otimes {\mathit{C}}^{\left(\mathrm{0}{g}_{\mathrm{2}}\mathrm{0}\right)}\otimes {\mathit{C}}^{\left(\mathrm{00}{g}_{\mathrm{3}}\right)}\right]\mathit{u}(15)

其中：\otimes表示克罗内克积运算符；{\mathit{C}}^{\left(g_{\mathrm{1}}\mathrm{00}\right)}\mathrm{、}{\mathit{C}}^{\left(\mathrm{0}{g}_{\mathrm{2}}\mathrm{0}\right)}\mathrm{、}

{\mathit{C}}^{\left(\mathrm{00}{g}_{\mathrm{3}}\right)}分别为关于x，y，t的微分矩阵。

结合式(1)和式(10)可以得到

{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^l}\left[\begin{array}{l}{L}_i\left(x\right)M_j\left(y\right)T_k^{″}\left(t\right)-\\ c^{\mathrm{2}}\left(\begin{array}{l}{L}_i^{″}\left(x\right)M_j\left(y\right)T_k\left(t\right)\\ +L_i\left(x\right)M_j^{″}\left(y\right)T_k\left(t\right)\end{array}\right)\end{array}\right]u_{ijk}=f\left(x,y,t\right)(16)

令方程(16)在以下节点成立\left\{\left(x_r,y_s,t_p\right),r=\mathrm{1},\mathrm{2},\cdots ,m;s=\mathrm{1},\mathrm{2},\cdots ,n;p=\mathrm{1},\mathrm{2},\cdots ,l\right\}，那么

\begin{array}{l}{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^l}\left[\begin{array}{l}{L}_i\left(x_r\right)M_j\left(y_s\right)T_k^{″}\left(t_p\right)-\\ c^{\mathrm{2}}\left(\begin{array}{l}{L}_i^{″}\left(x_r\right)M_j\left(y_s\right)T_k\left(t_p\right)\\ +L_i\left(x_r\right)M_j^{″}\left(y_s\right)T_k\left(t_p\right)\end{array}\right)\end{array}\right]u_{ijk}=f\left(x_r,y_s,t_p\right)\\ \end{array}(17)

线性方程组(17)写成如下矩阵的形式

\left[\begin{array}{l}\left({\mathit{I}}_m\otimes {\mathit{I}}_n\right)\otimes {\mathit{C}}^{\left(\mathrm{002}\right)}-\\ c^{\mathrm{2}}\left({\mathit{C}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_n\right)\otimes {\mathit{I}}_l-\\ c^{\mathrm{2}}\left({\mathit{I}}_m\otimes {\mathit{C}}^{\left(\mathrm{020}\right)}\right)\otimes {\mathit{I}}_l\end{array}\right]\mathit{U}=\mathit{F}(18)

简记为

\mathit{L}\mathit{U}=\mathit{F}(19)

其中

\begin{array}{l}\mathit{L}=\left({\mathit{I}}_m\otimes {\mathit{I}}_n\right)\otimes {\mathit{C}}^{\left(\mathrm{002}\right)}-c^{\mathrm{2}}\left({\mathit{C}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_n\right)\\ \otimes {\mathit{I}}_l-c^{\mathrm{2}}\left({\mathit{I}}_m\otimes {\mathit{C}}^{\left(\mathrm{020}\right)}\right)\otimes {\mathit{I}}_l\end{array}

\mathit{U}={\left[\begin{array}{l}{u}_{\mathrm{111}},u_{\mathrm{112}},\cdots ,u_{\mathrm{11}l},u_{\mathrm{121}},u_{\mathrm{122}},\cdots ,u_{\mathrm{12}l},\cdots ,\\ \cdots ,u_{mn\mathrm{1}},u_{mn\mathrm{2}},\cdots ,u_{mnl}\end{array}\right]}^{\mathrm{T}}

\mathit{F}={\left[\begin{array}{l}{f}_{\mathrm{111}},f_{\mathrm{112}},\cdots ,f_{\mathrm{11}l},f_{\mathrm{121}},f_{\mathrm{122}},\cdots ,f_{\mathrm{12}l},\cdots ,\\ \cdots ,f_{mn\mathrm{1}},f_{mn\mathrm{2}},\cdots ,f_{mnl}\end{array}\right]}^{\mathrm{T}}

考虑二维波动方程(1)的初始条件

\begin{array}{l}u\left(x,y,\mathrm{0}\right)=\varphi \left(x,y\right)\\ u_t\left(x,y,\mathrm{0}\right)=\psi \left(x,y\right)\end{array}(20)

以及x的边界条件

\begin{array}{l}u\left(a,y,t\right)={\zeta }_{\mathrm{1}}\left(y,t\right)\\ u\left(b,y,t\right)={\zeta }_{\mathrm{2}}\left(y,t\right)\end{array}(21)

和y的边界条件

\begin{array}{l}u\left(x,c,t\right)={\chi }_{\mathrm{1}}\left(x,t\right)\\ u\left(x,d,t\right)={\chi }_{\mathrm{2}}\left(x,t\right)\end{array}(22)

利用重心有理插值式(10)来近似式(20)~(22)，可以得到初始条件和边界条件的离散形式，最终离散的矩阵形式如下：

\begin{array}{l}\left[{\mathit{I}}_m\otimes {\mathit{I}}_n\otimes {\mathit{I}}_{l\mathrm{1}}\right]\mathit{u}=\mathit{\varphi }\\ \left[{\mathit{I}}_m\otimes {\mathit{I}}_n\otimes {\mathit{C}}_{t\mathrm{1}}^{\left(\mathrm{1}\right)}\right]\mathit{u}=\mathit{\psi }\end{array}(23)

其中：

\begin{array}{l}\mathit{\varphi }=\left[\varphi \left(x_{\mathrm{1}},y_{\mathrm{1}}\right),\cdots ,\varphi \left(x_{\mathrm{1}},y_n\right),\cdots \right.{\left.,\cdots ,\varphi \left(x_m,y_n\right)\right]}^{\mathrm{T}}\\ \mathit{\psi }=\left[\psi \left(x_{\mathrm{1}},y_{\mathrm{1}}\right),\cdots ,\psi \left(x_{\mathrm{1}},y_n\right),\cdots \right.\left.,\cdots ,\psi \left(x_m,y_n\right)\right]\end{array}

{\mathit{I}}_{l\mathrm{1}}表示l阶单位矩阵的第1行，{\mathit{C}}_{t\mathrm{1}}^{\left(\mathrm{1}\right)}表示t方向上的一阶微分矩阵的第1行。

\begin{array}{l}\left[{\mathit{I}}_{m\mathrm{1}}\otimes {\mathit{I}}_n\otimes {\mathit{I}}_l\right]\mathit{u}={\mathit{\zeta }}_{\mathrm{1}}\left(y,t\right)\\ \left[{\mathit{I}}_{mm}\otimes {\mathit{I}}_n\otimes {\mathit{I}}_l\right]\mathit{u}={\mathit{\zeta }}_{\mathrm{2}}\left(y,t\right)\end{array}(24)

\begin{array}{l}\left[{\mathit{I}}_m\otimes {\mathit{I}}_{n\mathrm{1}}\otimes {\mathit{I}}_l\right]\mathit{u}={\chi }_{\mathrm{1}}\left(x,t\right)\\ \left[{\mathit{I}}_m\otimes {\mathit{I}}_{nn}\otimes {\mathit{I}}_l\right]\mathit{u}={\chi }_{\mathrm{2}}\left(x,t\right)\end{array}(25)

其中：

\begin{array}{l}\left\{{\mathit{\zeta }}_p={\left[\begin{array}{l}{\zeta }_p\left(y_{\mathrm{1}},t_{\mathrm{1}}\right),\cdots ,{\zeta }_p\left(y_{\mathrm{1}},t_l\right),\\ \cdots ,\cdots ,{\zeta }_p\left(y_n,t_l\right)\end{array}\right]}^{\mathrm{T}},p=\mathrm{1},\mathrm{2}\right\}\\ \left\{{\mathit{\chi }}_q={\left[\begin{array}{l}{\chi }_q\left(x_{\mathrm{1}},t_{\mathrm{1}}\right),\cdots ,{\chi }_q\left(x_{\mathrm{1}},t_l\right),\\ \cdots ,\cdots ,{\chi }_q\left(x_m,t_l\right)\end{array}\right]}^{\mathrm{T}},q=\mathrm{1},\mathrm{2}\right\}\end{array}{\mathit{I}}_{m\mathrm{1}}表示m阶单位矩阵的第1行，{\mathit{I}}_{mm}表示m阶单位矩阵的第m行，{\mathit{I}}_{n\mathrm{1}}表示n阶单位矩阵的第1行，{\mathit{I}}_{nn}表示n阶单位矩阵的第n行。

最终得到二维波动方程(1)的离散矩阵(19)以及初边界条件的离散形式(23)~(25)，另外，本文采用置换法来计算边界条件。

4 数值算例

利用给定解析解的数值算例进行数值模拟来验证本文所提方法的高精度和高效性。

所有算例所选插值节点为第二类Chebyshev节点，

\begin{array}{l}{x}_i=-\mathrm{c}\mathrm{o}\mathrm{s}\frac{i\mathrm{\pi }}{m},i=\mathrm{0},\mathrm{1},\cdots ,m\\ y_j=-\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\mathrm{\pi }}{n},j=\mathrm{0},\mathrm{1},\cdots ,n\\ t_k=-\mathrm{c}\mathrm{o}\mathrm{s}\frac{k\mathrm{\pi }}{l},k=\mathrm{0},\mathrm{1},\cdots ,l\end{array}

以及等距节点

\begin{array}{l}{x}_i=a+i\frac{b-a}{m},i=\mathrm{0},\mathrm{1},\cdots ,m\\ y_j=c+j\frac{d-c}{n},j=\mathrm{0},\mathrm{1},\cdots ,n\\ t_k=k\frac{T}{l},k=\mathrm{0},\mathrm{1},\cdots ,l\end{array}

定义绝对误差E_{\mathrm{a}}和相对误差E_{\mathrm{r}}：

E_{\mathrm{a}}=||u_{\mathrm{c}}-u_{\mathrm{e}}|{|}_{\mathrm{2}},E_{\mathrm{r}}=\frac{||u_{\mathrm{c}}-u_{\mathrm{e}}|{|}_{\mathrm{2}}}{||u_{\mathrm{e}}|{|}_{\mathrm{2}}}

其中：{‖·‖}_{\mathrm{2}}表示向量的\mathrm{2}-范数；u_{\mathrm{c}}表示所求得的数值解；u_{\mathrm{e}}表示方程已知的精确解。

算例1 考虑如下常系数二维波动方程

\left\{\begin{array}{l}{u}_{tt}=u_{xx}+u_{yy}+f(x,y,t),\left(x,y,t\right)\in \Omega \times \left.\left(\mathrm{0},T\right.\right]\\ \begin{array}{l}u\left(x,y,\mathrm{0}\right)=x(x-\mathrm{1})+y(y-\mathrm{1})\\ u_t\left(x,y,\mathrm{0}\right)=\mathrm{\pi }\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)\end{array}\\ \begin{array}{l}u\left(\mathrm{0},y,t\right)=y(y-\mathrm{1})+\mathrm{2}{t}^{\mathrm{2}}\\ u\left(\mathrm{1},y,t\right)=y(y-\mathrm{1})+\mathrm{2}{t}^{\mathrm{2}}\end{array}\\ \begin{array}{l}u\left(x,\mathrm{0},t\right)=x(x-\mathrm{1})+\mathrm{2}{t}^{\mathrm{2}}\\ u\left(x,\mathrm{1},t\right)=x(x-\mathrm{1})+\mathrm{2}{t}^{\mathrm{2}}\end{array}\end{array}\right.

其中：\Omega =\left\{\left(x,y\right)|\mathrm{0}\le x,y\le \mathrm{1}\right\}，f(x,y,t)={\mathrm{\pi }}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)

\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }t)。本算例的精确解为

\begin{array}{l}u\left(x,y,t\right)=x(x-\mathrm{1})+y(y-\mathrm{1})+\mathrm{2}{t}^{\mathrm{2}}\\ +\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }t)\end{array}

取T=\mathrm{1}，在整个计算区域内分别选取\mathrm{16}\times \mathrm{16}\times \mathrm{16}个张量积型第二类Chebyshev插值节点以及等距插值节点来离散算例1中的波动方程，采用重心有理插值法和重心Lagrange插值法计算，将两种插值配点法所得到的结果进行比较。在使用重心有理插值法计算时，取参数d=\mathrm{10}。

图1给出了采用第二类Chebyshev节点在两种重心插值配点法下进行数值计算得到算例1的数值解和绝对误差。图1(a)、1(b)表明Chebyshev节点下的两种方法得到的数值解和精确解具有较高的一致性。图1(c)、1(d)表明两种方法计算的数值精度相差不大，精度量级均保持在10^-10。

表1给出t=\mathrm{1}时算例1在不同插值节点下的两种重心插值配点法的误差。结果表明，对于不同的插值节点，本文所提方法得到的数值解的误差和重心Lagrange插值法得到的数值解的精度一致，也就是说，两种插值法的数值拟合效果相当。对于重心有理插值配点法，Chebyshev节点下的重心有理插值配点法的效果略优于等距节点。

表2给出不同插值节点数下的两种方法的绝对误差，从表2中得知，随着插值节点数的增加，重心有理插值法在两种插值节点下始终保持良好的数值精度，在节点数达到\mathrm{17}\times \mathrm{17}\times \mathrm{17}时，能够保持较稳定的精度；重心Lagrange插值在Chebyshev节点下能够保持稳定的数值精度，在节点数达到\mathrm{13}\times \mathrm{13}\times \mathrm{13}时，计算精度相对较高，量级最终保持在{\mathrm{10}}^{-\mathrm{8}}，但是在等距节点下，重心Lagrange插值随着插值节点数增加到\mathrm{15}\times \mathrm{15}\times \mathrm{15}时，精度却逐渐下降，且下降幅度较大，在节点数为\mathrm{19}\times \mathrm{19}\times \mathrm{19}时，数值模拟极不稳定，重心有理插值法比重心Lagrange插值法在等距节点下稳定。

算例2 考虑如下二维波动方程

\left\{\begin{array}{l}{u}_{tt}=u_{xx}+u_{yy},\left(x,y,t\right)\in \Omega \times \left.\left(\mathrm{0},T\right.\right]\\ \begin{array}{l}u\left(x,y,\mathrm{0}\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(-x-y\right)\\ u_t\left(x,y,\mathrm{0}\right)=-\sqrt[]{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(-x-y\right)\end{array}\\ \begin{array}{l}u\left(\mathrm{0},y,t\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\sqrt[]{\mathrm{2}}t-y\right)\\ u\left(\mathrm{1},y,t\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\sqrt[]{\mathrm{2}}t-\mathrm{1}-y\right)\end{array}\\ \begin{array}{l}u\left(x,\mathrm{0},t\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\sqrt[]{\mathrm{2}}t-x\right)\\ u\left(x,\mathrm{1},t\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\sqrt[]{\mathrm{2}}t-\mathrm{1}-x\right)\end{array}\end{array}\right.

其中：\Omega =\left\{\left(x,y\right)|\mathrm{0}\le x,y\le \mathrm{1}\right\}，T=\mathrm{1}。本算例的精确解为

u\left(x,y,t\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\sqrt[]{\mathrm{2}}t-x-y\right)\mathrm{。}

采用\mathrm{14}\times \mathrm{14}\times \mathrm{14}个第二类Chebyshev节点对算例2中的方程进行离散，取d=\mathrm{10}。分别使用两种重心插值配点法进行数值计算，得到图2中的数值解和误差分布。图2(a)、2(b)表明，两种方法所得的数值解与精确解拟合度都较好。图2(c)、2(d)给出两种方法的绝对误差，表明在Chebyshev节点下，两种方法的误差精度量级均保持在{\mathrm{10}}^{-\mathrm{11}}。

表3给出算例2在不同插值节点下的两种重心插值配点法的绝对误差和相对误差，结果表明，采用Chebyshev节点的有理插值法的数值精度稍微高。选取等距节点时，重心有理插值法比重心Lagrange插值法拟合效果更好。

表4给出算例2在不同插值节点数下两种方法的绝对误差，从表4中得知，两种方法随着插值节点数的增加，数值计算精度虽然都有所下降，但是在Chebyshev节点下都能保持良好的数值精度。等距节点下的重心有理插值法计算精度虽有所下降，但是仍能保持良好的稳定性。对于重心Lagrange插值法，等距节点下该方法随着节点数增加到\mathrm{17}\times \mathrm{17}\times \mathrm{17}时，计算精度下降幅度大，达到了{\mathrm{10}}^{-\mathrm{1}}，节点数增加到\mathrm{19}\times \mathrm{19}\times \mathrm{19}时，该方法的误差变得非常大，说明此方法在等距节点下极不稳定。

算例3 考虑如下常系数二维波动方程

\left\{\begin{array}{l}{u}_{tt}=u_{xx}+u_{yy}+f(x,y,t),\left(x,y,t\right)\in \Omega \times \left.\left(\mathrm{0},T\right.\right]\\ \begin{array}{l}u\left(x,y,\mathrm{0}\right)=\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)\\ u_t\left(x,y,\mathrm{0}\right)=-\mathrm{\pi }\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)\end{array}\\ u\left(\mathrm{0},y,t\right)=u\left(\mathrm{1},y,t\right)=\mathrm{0}\\ u\left(x,\mathrm{0},t\right)=u\left(x,\mathrm{1},t\right)=\mathrm{0}\end{array}\right.

其中：\Omega =\left\{\left(x,y\right)|\mathrm{0}\le x,y\le \mathrm{1}\right\}，f(x,y,t)=\mathrm{3}{\mathrm{\pi }}^{\mathrm{2}}{\mathrm{e}}^{-\mathrm{\pi }t}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)

\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)。本算例的精确解为

u\left(x,y,t\right)={\mathrm{e}}^{-\mathrm{\pi }t}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }x)\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }y)。

取T=\mathrm{1}，采用\mathrm{20}\times \mathrm{20}\times \mathrm{20}个第二类Chebyshev节点对算例3中的方程进行离散，同样地，取d=\mathrm{10}。分别使用两种重心插值配点法进行数值计算，得到图3中的数值解和误差分布。图3(a)、3(b)表明，两种方法所得的数值解与精确解拟合度都较好。图3(c)、3(d)给出两种方法的绝对误差，表明在Chebyshev节点下两种方法的误差精度量级均保持在{\mathrm{10}}^{-\mathrm{11}}。

表5给出算例3在不同插值节点数下的两种重心插值配点法的绝对误差和相对误差，结果表明，采用Chebyshev节点的重心Lagrange插值法的数值精度稍高。选取等距节点时，重心有理插值法比重心Lagrange插值法精度好。

表6给出算例3两种方法在不同插值节点数下的绝对误差，从表6看出：两种方法随着插值节点数的增加，数值计算精度都有所上升，重心Lagrange插值在Chebyshev节点下的最高精度达到{\mathrm{10}}^{-\mathrm{10}}，重心有理插值法计算精度略低，为{\mathrm{10}}^{-\mathrm{9}}。等距节点下的重心有理插值法能够保持较为稳定的数值结果，但是对于重心Lagrange插值法，等距节点下该方法随着节点数增加到\mathrm{17}\times \mathrm{17}\times \mathrm{17}时，计算精度下降，达到{\mathrm{10}}^{-\mathrm{2}}，节点数增加到\mathrm{19}\times \mathrm{19}\times \mathrm{19}时，该方法的误差变得非常大，此方法在等距节点下不稳定。

5 结语

本文采用重心有理插值配点法求解二维波动方程，并且和重心Lagrange插值配点法进行了比较。数值算例表明，两种方法在Chebyshev节点下的数值解和精确解均具有较好的吻合度。在等距节点下，本文介绍的方法计算精度要高一些。

参考文献

[1] SCHOEDER S, KRONBICHLER M, WALL W A. Arbitrary high-order explicit hybridizable discontinuous Galerkin methods for the acoustic wave equation[J]. Journal of Scientific Computing, 2018, 76(2): 969-1006.

[2] ANTONIETTI P F, MAZZIERI I, MUHR M, et al. A high-order discontinuous Galerkin method for nonlinear sound waves[J]. Journal of Computational Physics, 2020, 415: 109484.

[3] YUE C, HIGAZY M, KHATER O M A, et al. Computational and numerical simulations of the wave propagation in nonlinear media with dispersion processes[J]. AIP Advances, 2023, 13(3): 035232.

[4] KUMAR D, KUMAR S. Some new periodic solitary wave solutions of (3+1)-dimensional generalized shallow water wave equation by lie symmetry approach[J]. Computers & Mathematics with Applications, 2019, 78(3): 857-877.

[5] ZHANG M, YAN W J, JING F F, et al. Discontinuous Galerkin method for the diffusive-viscous wave equation[J]. Applied Numerical Mathematics, 2023, 183: 118-139.

[6] YI N, HUANG Y, LIU H. A conservative discontinuous Galerkin method for nonlinear electromagnetic Schrödinger equations[J]. SIAM Journal on Scientific Computing, 2019, 41(6): B1389-B1411.

[7] HARARI I, TURKEL E. Accurate finite difference methods for time-harmonic wave propagation[J]. Journal of Computational Physics, 1995, 119(2): 252-270.

[8] ACHOURI T. Finite difference schemes for the two-dimensional semilinear wave equation[J]. Numerical Methods for Partial Differential Equations, 2019, 35(1): 200-221.

[9] WANG X P, GAO F Z, SUN Z J. Weak Galerkin finite element method for viscoelastic wave equations[J]. Journal of Computational and Applied Mathematics, 2020, 375: 112816.

[10] GHOUDI T, MOHAMED M S, SEAID M. Novel adaptive finite volume method on unstructured meshes for time-domain wave scattering and diffraction[J]. Computers & Mathematics with Applications, 2023, 141: 54-66.

[11] KHASI M, RASHIDINIA J, RASOULIZADEH M N. Fast computing approaches based on a bilinear pseudo-spectral method for nonlinear acoustic wave equations[J]. SIAM Journal on Scientific Computing, 2023, 45(4): B413-B439.

[12] CUI M, JI C C, DAI W. A finite difference method for solving the wave equation with fractional damping[J]. Mathematical and Computational Applications, 2023, 29(1): 2.

[13] HAN W M, SONG C H, WANG F, et al. Numerical analysis of the diffusive-viscous wave equation[J]. Computers & Mathematics With Applications, 2021, 102: 54-64.

[14] 王兆清, 李淑萍, 唐炳涛. 一维重心型插值: 公式、算法和应用[J]. 山东建筑大学学报, 2007, 22(5): 448-453.

[15] 马燕, 王兆清, 唐炳涛. 波动问题的高精度重心有理插值配点法[J]. 山东科学, 2012, 25(3): 80-87.

[16] 宋灵宇, 武莉莉, 卢梦双. 重心插值配点法求解二维Sobolev方程[J]. 西北师范大学学报(自然科学版), 2021, 57(6): 31-37, 44.

[17] 屈金铮, 李金, 苏晓宁. 求解非线性伪抛物方程的重心Lagrange插值配点法[J]. 山东大学学报(理学版), 2023, 58(4): 29-39.

[18] 袁洪旺, 王希胤, 李金. 波动方程的高精度数值解方法[J]. 计算物理, 2024, 41(4): 426-439.

[19] ZHAO Q L, WANG R W, WANG Z Q, et al. Barycentric rational collocation method for burgers' equation[J]. Journal of Mathematics, 2022, 2022(1): 2177998.

[20] LI J, QU J Z. Barycentric lagrange interpolation collocation method for solving the sine-gordon equation[J]. Wave Motion, 2023, 120: 103159.

基金资助

国家自然科学青年基金资助项目(12101482)

陕西省科技厅重点研发一般资助项目(2024M722604)