線性Poisson-Boltzmann方程基于L2投影的虛單元方法

陳鍵鏵， 陽 鶯

(桂林電子科技大學 數學與計算科學學院，廣西 桂林 541004)

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

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

針對線性 PBE，利用L2投影算子，設計了一種虛單元離散格式，給出了在多邊形網格上的虛單元求解，并給出H1誤差分析。數值實驗結果驗證了理論結果的正確性。

1 預備知識

考慮如下一類線性PBE：

(1)

其中：

Ωm為分子區域，Ωs為溶劑區域，εs=80，εm=2，k2為Debye-Huckel參數。

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

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

(2)

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

(3)

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

(4)

其中：

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

(5)

(6)

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

(7)

2 虛單元空間

對任意h和E∈Γh，假設存在一個正常數λ，滿足[6,10]：

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

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

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

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

(8)

對應的全局空間為

(9)

2個投影算子的定義為

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

(10)

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

?p1∈P1(E)。

(11)

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

并滿足如下連續性和正定性：

?u,v∈H1(Ω)。

(12)

在單元E∈Γh上的內積形式定義[11]為

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

(13)

(14)

(15)

其中，對于SE，假設存在2個正常數α*和α*，滿足

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

對于任意的uh,vh∈Vh, 定義

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

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

(16)

引理1[11]對任意E∈Γh和單元E上任意函數φ，存在一個正常數C，滿足

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

引理2[11]對任意的v∈Vh，任意整數k≥1，有

ah(v,v)≥

(17)

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

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

(18)

其中常數C不依賴于h。

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

?E∈Γh。

(19)

3 誤差估計

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

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

(20)

其中C為與h無關的常數。

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

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

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

(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)，有以下估計：

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

I1+I2+I3,

(25)

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

(26)

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

(27)

對于I3估計，有

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 數值實驗

用數值例子驗證虛單元法的有效性，基于文獻[12-13]代碼進行計算。

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

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

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

圖1 三角形、四邊形及五邊形組成的混合多邊形網格

圖2 L2模誤差

圖3 H1模誤差

5 結束語

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