陳鍵鏵， 陽(yáng) 鶯

(桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004)

Poisson-Boltzmann equation(PBE)是用來(lái)描述電荷密度分布和離子濃度的偏微分方程，廣泛應(yīng)用于物理、生物和化學(xué)等[1]領(lǐng)域。PBE是一類帶有分片常數(shù)和奇性的問題，常用的求解方法有有限差分法[2-3]和有限元法[4-5]等。在實(shí)際問題中，有限差分法和有限元法對(duì)于求解不規(guī)則界面時(shí)效果不佳。虛單元法由于使用多邊形網(wǎng)格，更適合求解不規(guī)則界面問題。

近年來(lái)，虛單元方法被廣泛應(yīng)用于各種偏微分方程的數(shù)值求解，最初在2012年由Brezzi等[6]提出，起初被應(yīng)用于Poisson方程，后來(lái)被推廣到其他方程[7-8]。虛單元方法對(duì)于網(wǎng)格方面的要求較低，可以比較靈活地應(yīng)用于不規(guī)則的多邊形網(wǎng)格，包括非凸多邊形網(wǎng)格。虛單元法的核心是虛單元空間中試探函數(shù)和檢驗(yàn)函數(shù)無(wú)顯式表達(dá)式，事實(shí)上，虛單元法只需局部離散函數(shù)空間的多項(xiàng)式子空間的知識(shí)來(lái)提供穩(wěn)定和精確的數(shù)值方法。通過(guò)引入合適的投影算子，將虛單元空間的函數(shù)及其梯度映射到多項(xiàng)式上，從而實(shí)現(xiàn)投影算子僅使用虛單元空間的自由度即可計(jì)算，而剛度矩陣便是通過(guò)自由度來(lái)進(jìn)行計(jì)算的。

針對(duì)線性 PBE，利用L2投影算子，設(shè)計(jì)了一種虛單元離散格式，給出了在多邊形網(wǎng)格上的虛單元求解，并給出H1誤差分析。數(shù)值實(shí)驗(yàn)結(jié)果驗(yàn)證了理論結(jié)果的正確性。

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

考慮如下一類線性PBE：

(1)

其中：

Ωm為分子區(qū)域，Ωs為溶劑區(qū)域，εs=80，εm=2，k2為Debye-Huckel參數(shù)。

其中：kB為Boltzmann常數(shù)；T為溫度；qi為第i中離子的電荷量；ec為單位電荷。δi(x)=δ(x-xi)為xi處的Dirac分布。

方程(1)具有奇性，為了避免求解奇性問題， 根據(jù)文獻(xiàn)[4]作如下分解：

(2)

將式(2)代入式(1)，得到以下PBE：

(3)

方程(3)的弱解為u∈H1(Ω)，滿足

(4)

其中：

a(u,v)=(εu,v)；

(5)

(6)

(f,v)=(·((ε-εm)G),v)。

(7)

2 虛單元空間

對(duì)任意h和E∈Γh，假設(shè)存在一個(gè)正常數(shù)λ，滿足[6,10]：

1)單元E相對(duì)于半徑為λhE的球是星型；

2)單元E每條邊的長(zhǎng)度大于或等于λhE。

考慮一階的虛單元空間，首先定義E上的局部空間。對(duì)?E∈Γh，

Δv|E∈P-1(E)},P-1(E)={0},

(8)

對(duì)應(yīng)的全局空間為

(9)

2個(gè)投影算子的定義為

(Π0vh-vh，p1)E=0,?p1∈P1(E)；

(10)

((Πvh)-vh,p1)E=0,

?p1∈P1(E)。

(11)

雙線性形式a(u,v)在單元E上的表示為

并滿足如下連續(xù)性和正定性：

?u,v∈H1(Ω)。

(12)

在單元E∈Γh上的內(nèi)積形式定義[11]為

SE((I-Π)uh,(I-Π)vh),

(13)

(14)

(15)

其中，對(duì)于SE，假設(shè)存在2個(gè)正常數(shù)α*和α*，滿足

α*aE(v,v)≤SE(v,v)≤α*aE(v,v),

對(duì)于任意的uh,vh∈Vh, 定義

方程(3)的虛單元法離散格式：存在uh∈Vh，使得

ah(uh,vh)+bG,h(uh,vh)=(fh,vh),?vh∈Vh。

(16)

引理1[11]對(duì)任意E∈Γh和單元E上任意函數(shù)φ，存在一個(gè)正常數(shù)C，滿足

m、s∈N,m≤s≤k+1,s≥1；|φ-φI‖m,E≤

引理2[11]對(duì)任意的v∈Vh，任意整數(shù)k≥1，有

ah(v,v)≥

(17)

引理3[11]雙線性形式ah(·,·)在Vh×Vh上是連續(xù)的，即

ah(uh,vh)≤C‖uh‖1‖vh‖1,uh,vh∈Vh,

(18)

其中常數(shù)C不依賴于h。

引理4[11]對(duì)任意的u∈H2(Ω)和vh∈Vh， 有

?E∈Γh。

(19)

3 誤差估計(jì)

定理1假設(shè)u∈H2(E)∩W1,∞(E)是式(4)的解，f∈H1(E)，uh∈Vh是式(16)的解，則有

‖u-uh‖1≤C(h+‖u-uh‖0),

(20)

其中C為與h無(wú)關(guān)的常數(shù)。

證明將誤差e拆寫成2部分，

e=u-uh=u-uI+uI-uh,

記eh=uh-uI。對(duì)式(17)運(yùn)用龐加萊不等式，有

(fh,eh)-bG,h(uh,eh)-ah(uI,eh),

(21)

-ah(uI,eh)=ah(Π0u-uI,eh)-ah(Π0u,eh)=

ah(Π0u-uI,eh)-ah(Π0u,eh)+

a(Π0u,eh)+a(u-Π0u,eh)-a(u,eh)。

(22)

將式(22)中的a(u,eh)用式(4)替換，再將式(22)代入式(21)，得

bG,h(uh,eh))+(a(Π0u,eh)-ah(Π0u,eh))+

ah(Π0u-uI,eh)+a(u-Π0u,eh)=

H1+H2+H3+H4+H5。

(23)

由式(15)，有以下近似：

(24)

由式(14)，有以下估計(jì)：

H2=bG(u,eh)-bG,h(uh,eh)=

I1+I2+I3,

(25)

C‖u-Π0u‖0‖eh‖0≤Ch‖u‖1‖eh‖。

(26)

由文獻(xiàn)[3]，u∈L∞(Ω),G∈C∞(Ωs)，有

(27)

對(duì)于I3估計(jì)，有

C‖Π0u-Π0uh‖0‖Π0eh‖0≤

C‖u-uh‖0‖eh‖1。

(28)

將式(26)～(28)代入式(25)，可得

H2≤C(h(‖u‖1+1)+‖u-uh‖0)‖eh‖1≤

C(h+‖u-uh‖0)‖eh‖1。

(29)

由式(19)可得

H3=a(Π0u,eh)-ah(Π0u,eh)=

(30)

由式(18)得

H4=ah(Π0u-uI,eh)≤

‖Π0u-uI‖1‖eh‖1≤Ch‖eh‖1。

(31)

由式(12)，有

H5=a(u-Π0u,eh)≤

‖u-Π0u‖1‖eh‖1≤Ch‖eh‖1。

(32)

將式(24)、(29)～(32)代入式(23)， 可推得式(20)成立。證畢。

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

用數(shù)值例子驗(yàn)證虛單元法的有效性，基于文獻(xiàn)[12-13]代碼進(jìn)行計(jì)算。

圖1為三角形四邊形五邊形混合的多邊形網(wǎng)格。圖 2、3為虛單元法在圖1網(wǎng)格上的解與數(shù)值解的L2和H1模誤差圖像，該圖像在對(duì)數(shù)尺度下給出。根據(jù)文獻(xiàn)[14]，分別給出L2、H1模定義：

‖Πu*-Πuh‖0，|Πu*-Πuh|1。

從圖2、3可看出，L2模誤差達(dá)到二階，H1模誤差達(dá)到一階，與理論相符。由于問題(4)的解很小，圖2、3給出的數(shù)據(jù)是原問題的1018倍。

圖1 三角形、四邊形及五邊形組成的混合多邊形網(wǎng)格

圖2 L2模誤差

圖3 H1模誤差

5 結(jié)束語(yǔ)

針對(duì)線性PBE，利用L2投影算子，設(shè)計(jì)了一種虛單元離散格式，給出了H1范數(shù)的誤差分析。數(shù)值實(shí)驗(yàn)結(jié)果表明了虛單元法的有效性。這一方法可進(jìn)一步推廣到求解非線性PBE問題中。