泊松方程的有限差分方法及快速實現(xiàn)

劉 昊，張榮培，霍俊蓉

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

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

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

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

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

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

1 預(yù)備知識

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

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

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

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

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

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

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

定理2對矩陣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）中的拉普拉斯算子進行泰勒展開：

離散后為

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

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

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

令K1d=-Gx，則

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

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

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

上式可由3步來實現(xiàn)：

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

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

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

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

2.2 二維情形

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

式（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)；邊界項W=vec(Ux1+Uy1)。

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

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

對Gy?Ix可做類似運算，則

上式可由以下3步來實現(xiàn)：

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

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

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

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

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

2.3 三維情形

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

式（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由邊界值來確定，且W=vec(Ux1+Uy1+Uz1)。

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

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

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

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

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

由式（18）分解可得

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

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

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

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

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

3 數(shù)值實驗

針對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 一維情形下使用快速離散正弦變換和使用追趕法的誤差、誤差階及時間比較

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

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

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

4 結(jié)語

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