二维泊松方程基于离散正弦变换的高阶紧致差分方法

邹志涵，向远强

（贵州师范大学 数学科学学院，贵州 贵阳 550025）

椭圆型泊松方程在电磁学、理论物理和工程计算方面有着广泛的应用，许多典型问题都可以归纳为泊松方程边值问题。由于计算的复杂性，泊松方程的公式解并不易得到。随着计算机软硬件技术的发展，计算其数值解变得越来越简单。求解偏微分方程通常有3种方法，有限差分法是最基本的求解方法，其实质就是运用离散化方法求解微分方程，但有限差分格式的精度通常只有2阶。因此，开发高阶差分格式的有限差分法对于精确求解偏微分方程具有重要的意义。

与传统的差分格式相比，高阶紧致差分格式有模板短、精度高等优点。自从KREISS等［1］提出高阶紧致差分格式以来，人们已开发出多种紧致差分格式，并将这些紧致差分格式用于求解各类偏微分方程边值问题，得到了高精度数值解［2–9］。紧致差分格式的优点是使用少量的点就可以获得高精度解，缺点是求解时产生了矩阵系统。为了避免求解矩阵系统，MOHEBBI等［10］利用4阶紧致差分格式和c样条函数求解了一维对流扩散方程。ZHANG等［11］将多对角矩阵分裂成多个三对角矩阵，降低了计算难度。LI等［12］将高阶紧致差分方法与4阶Runge–Kutta方法结合，提高了KdV-Burgers方程的数值解精度。以此为基础，我们在求解泊松方程的Dirichlet边界问题中，提出了基于快速离散正弦变换的8阶精度的紧致差分格式，利用8阶紧致差分格式对二维泊松方程进行了离散处理，并利用所提出的快速正弦离散方法求解离散后的线性系统，通过数值算例证明该算法的准确性和有效性。

1 8阶紧致差分格式

对于方程

其中α、β、a、b和c都是待确定的常数。LELE［13］提出了一种关于二阶导数的近似方法，即将式（1）中每一项关于u的函数进行泰勒展开，得到

将式（2）代入式（1），合并合同的导数项后，可得关于α、β、a、b和c的方程组

LELE认为，在式（3）的5个等式中按顺序提取部分或全部方程组可以构造不同精度的公式。例如：若满足式（3）第一个等式，就可以构造2阶精度的紧凑差分近似格式；若满足式（3）前两个等式，就可以构造4阶精度的紧凑差分近似格式。以此类推，可以构造10阶精度的紧凑差分近似格式。

在本文中，我们用前4个方程构造紧致差分格式。为了缩短格式的模板长度，令c=0,则当时，可得8阶紧致差分格式。

将α、β、a、b和c的值代入式（1）后，经过化简，可得

2 传统有限差分方法的求解步骤

现在以二维泊松方程为例给出传统有限差分方法的求解步骤。需要求解的二维泊松方程为

对二维泊松方程应用8阶紧致差分格式进行离散，经过整理可得

若用矩阵表示，有

其中， ⊗表示矩阵的Kronecker乘积，I1∈R(N−1)×(N−1),I2∈R(M−1)×(M−1)，A为8阶紧致差分格式右端项的系数所组成的 (M− 1)×(M−1)维的五对角矩阵。同理，B是由左端项系数组成的 (N− 1)×(N−1)维的五对角矩阵。

相应的矩阵形式为

其中，

由方程的边界条件可知u0,j=uM,j= 0,ui,0=ui,N=0,故可做以下假设：当 0≤j≤N时，= 0;0≤i≤M时，

联立式（7）和式（9），可得

用传统的有限差分方法求解偏微分方程的数值解U，需要求矩阵的逆矩阵，而求逆矩阵的过程是非常耗时的。随着M和N的增大，需要求解的线性方程组之规模也会增大，计算成本和储存空间将会大幅增加。因此，有效的计算方法既要避免求逆矩阵，又要降低计算成本。用快速离散正弦变换法构造相应的线性方程组，就能避免求逆矩阵，降低计算成本，节省存储空间。

3 快速离散正弦变换方法

在二维空间中，u的离散正弦变换格式及的逆变换的格式的定义为

其中：k=1,2, …,M− 1;l= 1,2, … ,N−1。

利用式（10）也可以近似得出式（4）和式（5）中所有关于函数的项。将式（4）中所有关于函数的项用离散正弦变换公式替代，并用三角函数的有关公式化简，可得

同理，将式中所有关于函数u的项用离散正弦变换公式替代，并用三角函数的有关公式化简，可得

二维泊松方程在 (xi,yj)处可以离散为式（8），对离散后的泊松方程作离散正弦变换，可得

其中：i= 1,2, …,M− 1;j= 1,2, … ,N−1。

将式（12）和式（13）代入式（14），可得

得到，而fi,j=f(xi,yj)是已知函数在离散点处的值。

利用快速离散正弦变换方法求二维泊松方程的步骤表示如下：

步骤1：给定M和N，用M和N剖分网格，将已知函数f(x,y) 代入网格点 (xi,yj)，得到F。

利用8阶紧致差分格式将泊松方程离散后，利用求逆矩阵求数值解的计算成本为O(M2N2)，利用快速离散正弦方法求数值解的计算成本为O(MNlog2MN)。由此可见，M和N越大，快速离散正弦方法的优势就越明显。另外，快速离散正弦变换方法中没有求逆矩阵这一步，占用内存更少，减少了计算存储量，缩短了计算时间。

4 数值实验

现在利用基于8阶紧致差分格式的快速离散正弦变换方法求解二维的泊松方程。

算例是带有Dirichlet边界条件的泊松方程边值问题，计算和测试过程均在Matlab2020a中实现，实验结果以表格形式给出。其中，相对误差为收敛阶为u表示数值解，uexact表示精确解。

例1：分析泊松方程的零Dirichlet边值问题

其精确解uexact=sin(x)sin(y)。

表1给出了不同网格步长下分别利用本文给出的8阶紧致差分格式和文献［8］给出的差分格式计算所得到的结果。从表1中可以看出：文献［8］和本文给出的差分格式都能达到理论上的精度；利用本文给出的紧致差分格式能得到8阶精度的结果；本文给出的格式在不同步长下的相对误差均小于文献［8］中格式的相对误差；随着网格数的增加，本文给出的差分格式的相对误差明显小于文献［8］给出的差分格式的相对误差。这说明用本文给出的方法求出的数值解与精确解的吻合程度比文献［8］的高。

表1 不同网格下例1的相对误差和收敛阶

表2给出了利用矩阵求逆方法和本文给出的快速离散正弦方法求解例1所用的时间。从表2可以看出：网格点越密集，快速离散正弦算法的计算时间与矩阵运算的计算时间之间的差距就越大，故在层面上，快速离散正弦变换算法在计算时间上具有明显的优势；当M=N=256时，矩阵求逆方法因内存不足而无法进行，但离散正弦方法的计算时间却保持在0.1 s以内。

表2 两种方法在计算时间上的对比

注：1.N.A表示计算机内存溢出无法进行计算；2.计算时间为5次计算时间的平均值。

例2：分析泊松方程非零Dirichlet边值问题

其精确解

令U(x,y) =u(x,y)−x−y,将非零边界条件转换为零边界条件，可得

运用本文给出的方法求转换后的零边值问题的近似解U(x,y) ，再利用u(x,y) =U(x,y) +x+y即可求出u(x,y)。

表3给出了不同网格步长下利用本文给出的8阶紧致差分格式与文献［8］给出的差分格式所得到的计算结果。从表3中可以看出：利用紧致差分格式求解泊松方程非零Dirichlet边值问题时，其收敛阶可以达到8阶精度，与理论精度基本一致；利用本文给出的差分格式进行计算，得到的相对误差小于利用文献［8］给出的格式得到的相对误差。

表3 不同网格下例2的相对误差和收敛阶

4 结束语

针对泊松方程的Dirichlet边值问题，我们提出了一种基于快速离散正弦变换的8阶精度的紧致差分格式。借助8阶精度的紧致差分格式对泊松方程进行离散，利用快速离散正弦方法求解离散后的线性系统，避免了矩阵运算，降低了计算难度和内存损耗。与求逆矩阵方法相比，本文给出的方法在计算时间上有着显著的优势。通过数值实验证明了本文给出的方法的可行性。本文给出的方法还可以用于求解高维泊松方程和更复杂的方程。