劉 昊，張榮培，霍俊蓉

（沈陽(yáng)師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，遼寧 沈陽(yáng) 110136）

泊松方程是常見(jiàn)于靜電學(xué)、機(jī)械工程和理論物理的偏微分方程，在無(wú)引力源的情況下，ΔΦ=0（即拉普拉斯方程）；當(dāng)考慮引力場(chǎng)時(shí)，ΔΦ=f（f為引力場(chǎng)的質(zhì)量分布）。該方程通常用格林函數(shù)法求解，也可以采用分離變量法和特征線法求解［1］。

帶有規(guī)則區(qū)域Dirichlet邊界條件的泊松方程：

式中，Δ 表示拉普拉斯算子；x∈Ω；Ω ?Rd；當(dāng)d=1時(shí)，Δu=uxx；當(dāng)d=2 時(shí)，Δu=uxx+uyy；當(dāng)d=3 時(shí)，Δu=uxx+uyy+uzz。

由于帶有Dirichlet 邊界的泊松方程的解析解不容易求出，一般情況下采用有限差分法、有限元法和有限體積方法的數(shù)值方法進(jìn)行求解［2-4］。但這些方法在求解過(guò)程中需要求解大型的非線性代數(shù)方程組，計(jì)算過(guò)程復(fù)雜，耗費(fèi)時(shí)間較多。

綜上所述，本文提出了一種格式構(gòu)造簡(jiǎn)單、發(fā)展比較成熟的有限差分方法，可求解一維、二維和三維情形下的泊松方程，針對(duì)每一種情形分別采用直接求解方法和快速離散正弦變換方法進(jìn)行求解。

1 預(yù)備知識(shí)

將一維、二維和三維情形的求解區(qū)域、網(wǎng)格步長(zhǎng)和數(shù)值解矩陣統(tǒng)一列在表1中。

表1 d維泊松方程的求解區(qū)域、網(wǎng)格步長(zhǎng)和數(shù)值解矩陣

對(duì)式（1）做有限差分離散，所得有限差分矩陣Gp如下：

式中，p=(x,y,z)。該矩陣可以做如下特征分解：

定理1有限差分方程的微分矩陣Gx可進(jìn)行對(duì)角化［6］：

式 中，Λx=diag。對(duì)Gy和Gz也可得類似的結(jié)果。Sx是正弦矩陣，Sx與向量v 做乘積相當(dāng)于對(duì)v進(jìn)行離散正弦變換，可以結(jié)合傅里葉變換來(lái)實(shí)現(xiàn)，使得計(jì)算復(fù)雜度降低為O(Nxlog2Nx)。此時(shí)，Sx還是一個(gè)正交矩陣，即Sx=Sx-1。

由于二維和三維情形的有限差分矩陣需要用Kronecker積（?）得到，為了將其對(duì)角特征分解，需要用到下列Kronecker積的性質(zhì)：

定理2對(duì)矩陣A、B、C、D以及Kronecker 積（?）的性質(zhì)［5］：

1）(A?B)(C?D)=AC?BD；

2）(A?B)-1=A-1?B-1；

3）(BT?A)Vec(X)=Vec(AXB)；

4）(A?B)T=AT?BT。

2 有限差分

2.1 一維情形

采用二階中心差分得到差分格式，將式（1）中的拉普拉斯算子進(jìn)行泰勒展開：

離散后為

式中，fi為源函數(shù)；i∈Nx；fi=f(xi)的兩個(gè)邊界值為g(xb)和g(xe)，放在右端項(xiàng)中。

如果式（1）在x、y、z方向的Dirichlet邊界，那么在邊界處不會(huì)存在奇性。由于式（1）的邊界條件滿足以上假設(shè)，故不存在邊界交叉點(diǎn)處的奇性。

將式（3）改寫成矩陣形式：

令K1d=-Gx，則

式中，右端向量F的定義如下：

式中，i∈(2,3,…,Nx-1)。

由于Gx是三對(duì)角矩陣，故可使用追趕法［6］求解式（4）。此外，還可以采用離散正弦變換方法求解。由式（3）可得

上式可由3步來(lái)實(shí)現(xiàn)：

1）V1=Sx-1F可采用逆離散正弦變換實(shí)現(xiàn)；

2）對(duì)于V2=Λx-1?V1，由于Λx是對(duì)角矩陣，所以只有Nx次乘除法的計(jì)算過(guò)程較為復(fù)雜；

3）V=-SxV2采用離散正弦變換實(shí)現(xiàn)。

由于追趕法的計(jì)算復(fù)雜度為5Nx-4［7］，但離散正弦變換方法中的1）步和3）步的計(jì)算復(fù)雜度大約為O(Nxlog2Nx)［7］。由此可見(jiàn)，一維情形追趕法的計(jì)算方法具有優(yōu)勢(shì)。

2.2 二維情形

采用二階中心差分算子，將式（1）中的拉普拉斯算子進(jìn)行泰勒展開：

式（1）的差分格式如下

設(shè)Ip為Np階單位矩陣，將解矩陣U2向量化后得

根據(jù)Kronecker積的定義，式（5）的二階中心差分格式可以寫成如下形式：

式中，F(xiàn)=vec(F)；矩陣F的元素由源函數(shù)得到；Fi,j=f(xi,yj)；i∈(1,…,Nx)；j∈(1,…,Ny)；邊界項(xiàng)W=vec(Ux1+Uy1)。

Ux1和Uy1的元素由邊界值來(lái)確定：

式（8）同樣可以采用高斯消去法和離散正弦變換法求解。對(duì)K進(jìn)行特征分解，由于利用Kronecker積的性質(zhì)可得：

對(duì)Gy?Ix可做類似運(yùn)算，則

上式可由以下3步來(lái)實(shí)現(xiàn)：

1）利用Kronecker 積的性質(zhì)對(duì)式（8）右端項(xiàng)進(jìn)行拉直變換：

再應(yīng)用逆離散正弦變換實(shí)現(xiàn)：

2）對(duì)于V2=，./為矩陣中每個(gè)元素相除，這步計(jì)算的復(fù)雜度為Nx×Ny次乘除法；

3）U2=dst(dstV2)，即可直接采用離散正弦變換實(shí)現(xiàn)。

根據(jù)高斯消元法，當(dāng)Nx=Ny=N時(shí)，其計(jì)算復(fù)雜度為O(N4)，但快速離散正弦變換方法中1）步和3）步的計(jì)算復(fù)雜度僅為O(N3)［8-9］。由此可見(jiàn)，離散正弦變換方法的計(jì)算速度具有優(yōu)勢(shì)。

2.3 三維情形

三維情形下將式（1）中的拉普拉斯算子進(jìn)行泰勒展開：

式（1）的二階中心差分格式如下：

將解矩陣U3向量化后得

利用Kronecker 積定義，式（12）的差分格式可改寫為

式中，F(xiàn)i,j,k=f(xi,yj,zk)，i∈(1,…,Nx)；j∈(1,…,Ny)；k∈(1,…,Nz)；Ip為Np階單位矩陣；p=(x,y,z)；W由邊界值來(lái)確定，且W=vec(Ux1+Uy1+Uz1)。

Ux1、Uy1和Uz1的元素由邊界值確定：

故可將差分格式（14）改寫為

式（16）有兩種解法：第一種是利用高斯消元法直接求解式（15）；第二種利用離散正弦變換方法，對(duì)K3d進(jìn)行矩陣分解：

三維數(shù)組U3可以被視為二維數(shù)組上所有這樣的一維向量的集合，即

式中，U3(:,j,k)是通過(guò)固定最后兩個(gè)指標(biāo)j和k得到的Nx維向量。

由式（18）分解可得

式（21）可通過(guò)如下3步實(shí)現(xiàn)：

1）應(yīng)用逆離散正弦變換得到向量V1=，若將U3視為二維數(shù)組上的一維向量的集合，則有

3）直接應(yīng)用離散正弦變換實(shí)現(xiàn)，若將U3視為二維數(shù)組上的一維向量的集合，則有

由于積分的部分采用梯形求積公式進(jìn)行離散，故精度為二階精度。

在上述實(shí)現(xiàn)過(guò)程中，利用高斯消元法求解三維泊松方程時(shí)，其計(jì)算量為O(N7)［7］，但在快速離散正弦變換方法中，1）步和3）步的計(jì)算復(fù)雜度僅為O(N2log2N)［7］。由此可見(jiàn)，離散正弦變換方法的計(jì)算速度極具有優(yōu)勢(shì)。

3 數(shù)值實(shí)驗(yàn)

針對(duì)d維齊次泊松方程，一維情形下選取精確解函數(shù)為u(x)=π2sin(πx)，結(jié)果如表2 所示；二維情形下選取精確解函數(shù)為u(x,y)=2π2sin(πx)sin(πy)，結(jié)果如表3 所示；三維情形下選取精確解函數(shù)為u(x,y,z)=2π2sin(πx)sin(πy)sin(πz)，結(jié)果如表4所示。

表2 一維情形下使用快速離散正弦變換和使用追趕法的誤差、誤差階及時(shí)間比較

表3 二維情形下使用快速離散正弦變換和使用追趕法的誤差、誤差階及時(shí)間比較

表4 三維情形下使用快速離散正弦變換和使用追趕法的誤差、誤差階及時(shí)間比較

上述結(jié)果表明，本文提出的快速正弦離散變換方法具有更高的效率。

4 結(jié)語(yǔ)

本文應(yīng)用快速離散正弦變換求解泊松方程，該方法能夠在二維和三維的情況下快速地求解離散后得到的非線性代數(shù)方程組，提高了計(jì)算效率。實(shí)驗(yàn)結(jié)果表明：在一維情形時(shí)，在同樣的精確度下利用追趕法直接求解泊松方程的計(jì)算時(shí)間較快，正弦離散變換方法的速度較慢；但在二維和三維的情形下，快速正弦離散變換法的計(jì)算時(shí)間更短，計(jì)算次數(shù)更少，計(jì)算復(fù)雜度更低，且同樣可以達(dá)到較高的精確度。