孫高峰,牛瑞萍,李 明

(太原理工大學(xué) 數(shù)學(xué)學(xué)院,太原 030024)

復(fù)雜生物膜動(dòng)力系統(tǒng)的顯式有限差分方法研究

孫高峰,牛瑞萍,李 明

(太原理工大學(xué) 數(shù)學(xué)學(xué)院,太原 030024)

針對刻畫生物膜成長的動(dòng)力系統(tǒng)設(shè)計(jì)了一個(gè)非負(fù)、有界的顯式差分格式,系統(tǒng)是由四個(gè)變系數(shù)的非線性拋物型偏微分方程組成,描述了營養(yǎng)基質(zhì)上多種微生物成分間的交互過程；所設(shè)計(jì)的有限差分格式巧妙地將非線性項(xiàng)轉(zhuǎn)化成線性形式,確保得到的數(shù)值解是非負(fù)、有界的。最后通過兩個(gè)數(shù)值算例檢驗(yàn)了方法的性能,并給出了宏觀意義下的數(shù)值解的誤差。

復(fù)雜生物膜系統(tǒng);非線性反應(yīng)擴(kuò)散方程組;顯式有限差分法;非負(fù)、有界

在科學(xué)計(jì)算中經(jīng)常要用數(shù)值求解各類微分方程[1],有限元[2-3]、有限差分法[6-8]和有限體積法[4-5]是最重要的常用方法。本文使用有限差分法求解生物膜系統(tǒng)的數(shù)值解.生物膜系統(tǒng)本身物理屬性要求數(shù)值解是非負(fù)、有界的,因此,設(shè)計(jì)非負(fù)、有界的有限差分技術(shù)已成為計(jì)算數(shù)學(xué)的重要研究領(lǐng)域[6-8]。

Eberl et al建立的拋物型偏微分方程組來描述復(fù)雜的生物系統(tǒng)與周圍營養(yǎng)基的相互作并證明了該問題在合理的解析約束下,非負(fù)、有界解的存在唯一性[9-11]。受Eberl所做工作的啟發(fā),目前的工作主要是制定一種新的顯式有限差分法來有效地求解非線性拋物型偏微分方程組。由于生物膜系統(tǒng)是具有實(shí)際意義的體系,其未知變量表示各種微生物及營養(yǎng)成分的密度,是非負(fù)、有界的變量。如何設(shè)計(jì)一個(gè)有效的差分格式確保數(shù)值解的非負(fù)、有界性尤為重要。通常的生物系統(tǒng)的擴(kuò)散系數(shù)是給定非負(fù)常數(shù),這樣的系統(tǒng)設(shè)計(jì)數(shù)值計(jì)算相對容易,文獻(xiàn)中[6-11]已有許多成功的方法，但Eberl et al的生物系統(tǒng)的擴(kuò)散系數(shù)是未知變量的函數(shù),使得控制方程組是一個(gè)非線性的變系數(shù)偏微分方程組,這增加了設(shè)計(jì)數(shù)值計(jì)算方法的難度;本文的主要內(nèi)容是在初邊值條件下對復(fù)雜生物膜系統(tǒng)顯式有限差分方法進(jìn)行研究。

1 生物背景及數(shù)學(xué)模型

復(fù)雜生物膜系統(tǒng)由下面四個(gè)非線性偏微分方程確定:

(1)

(2)

(3)

(4)

M(X,t)=B(X,t)+I(X,t)+E(X,t) .

(5)

本文,假定μ,kL,kE,kI,kS,YE為非負(fù)常數(shù),YH正常數(shù),這些所有的常數(shù)均是通過實(shí)驗(yàn)獲得的。和·分別表示空間梯度算子和散度算子。S表示關(guān)于最大值標(biāo)準(zhǔn)化后營養(yǎng)濃度;函數(shù)M(X,t)表示問題域上的系統(tǒng)生物質(zhì)量密度。復(fù)雜生物膜系統(tǒng)中包含種類型的微生物:B表示活的微生物密度;I表示惰性或死的微生物密度;E表示胞外聚合物密度,這三個(gè)變量是關(guān)于最大量標(biāo)準(zhǔn)化的得到的量。顯然,對于和M∈[0,1)。從生物學(xué)角度看,偏微分方程組(1)-(5)刻畫的是不同生物之間以及生物和營養(yǎng)基之間的相互交互轉(zhuǎn)化過程。函數(shù)D(u)描述生物演化動(dòng)力學(xué)(生長、衰減、繁衍)規(guī)律的非線性擴(kuò)散系數(shù),通常是微生物的函數(shù),因此它是非線性的“根源”,本章使用下列形式的擴(kuò)散系數(shù):

式中,α≥1,β≥1.

方程(1)刻畫的是生物膜上的營養(yǎng)的消耗和轉(zhuǎn)化,值得注意的是活的微生物和胞外聚合物是營養(yǎng)的潛在“源”(可以轉(zhuǎn)化成營養(yǎng)),右端的第三和第四項(xiàng)刻畫的是轉(zhuǎn)化機(jī)制。這里,正常數(shù)kL表示活生物轉(zhuǎn)化為營養(yǎng)物的效率;正常數(shù)kE表示胞外聚合物轉(zhuǎn)化成營養(yǎng)物的效率;右端第一項(xiàng)營養(yǎng)物的擴(kuò)散項(xiàng);右端第二項(xiàng)是標(biāo)準(zhǔn)的Monod-type函數(shù),描述的是活生物如何消耗營養(yǎng)物.參數(shù)μ表示活生物的最大增長率;參數(shù)1/YH表示活生物消耗營養(yǎng)率。

方程(2)刻畫了生物膜上活生物的增長和轉(zhuǎn)。右端第一項(xiàng)控制活生物的擴(kuò)散;第二項(xiàng)是標(biāo)準(zhǔn)的Monod-type函數(shù),描述的是營養(yǎng)物如何促進(jìn)活生物的增長;第三和第四項(xiàng)描述了活生物如何轉(zhuǎn)化成營養(yǎng)物和惰性生物;正常數(shù)kI表示活生物轉(zhuǎn)化成惰性生物的效率。

方程(3)右端的第一項(xiàng)控制惰性生物的擴(kuò)散,第二項(xiàng)表示活生物如何促進(jìn)惰性生物的增長。

方程(4)描述的是胞外聚合物的消耗和轉(zhuǎn)化。方程右端的第一項(xiàng)控制胞外聚合物的擴(kuò)散;第二項(xiàng)描述的活生物如何促進(jìn)胞外聚合物的增長;常數(shù)YE表示胞外聚合物消耗活生物的效率。

2 有限差分方法

下面我們給出由方程式(1)-(5)所控制復(fù)雜生物膜系統(tǒng)的顯式非負(fù)、有界差分逼近形式。設(shè)問題域?yàn)檎?guī)的時(shí)-空區(qū)域:

(x,y,t)∈[0,X]×[0,Y]×[0,T] .

控制方程(1)-(4)的擴(kuò)散項(xiàng)可統(tǒng)一表示成下列方程:

·(D(M)f)=D(M)xfx+

D(M)yfy+D(M)(fxx+fyy) ,

式中：函數(shù)f依次可代表S,B,I,E。

下面我們使用以下有限差分形式對系統(tǒng)進(jìn)行近似:

δxxf=

δyyf=

為簡化記號,我們令

下面給出控制方程(1)-(4)式中擴(kuò)散項(xiàng)的逼近形式為:

·(D(M)f)≈

其中,函數(shù)f依次可代表S,B,I,E。為進(jìn)一步簡化記號,令

因此,復(fù)雜的生物膜系統(tǒng)(1)-(4)可由下面有限差分近似:

3 離散有限差分格式

為了表述的方便,下面給出一維模型的有限差分格式。

設(shè)a,b為實(shí)數(shù),且a

xi=(i-1)Δx.

i=1,2,3,…,M+1,且Δx=(b-a)/M.tk=(k-1)Δt.k=1,2,3,…,P+1,且Δt=T/P.

(6)

(7)

(8)

(9)

(10)

為了說明所提出的有限差分格式的優(yōu)越性,下面給出另外一種標(biāo)準(zhǔn)中心有限差分格式:

(11)

(12)

(13)

(14)

這里

4 數(shù)值算例

下面根據(jù)系統(tǒng)定義邊界為Neumann邊界條件的系統(tǒng)宏觀誤差。

聯(lián)立系統(tǒng)方程(1)-(5),有:

(15)

將上式方程在域(x,y,t)∈[0,X]×[0,Y]×[0,t]上積分,得到下列積分方程形式:

式中：b1=μc2,b2=(c1-c5)c2,b3=(c1-c5)(c2+c4)-μc2,b4=μc4.

進(jìn)一步簡化記號,令

這樣方程(15)可以改寫成下列形式:

b1(St-S0)+b2(Bt-B0)+

b3(It-I0)+b4(Et-E0)=0 .

因此，宏觀誤差Em可定義為下列形式:

Em=|b1(SNt-S0)+b2(BNt-B0)+

b3(INt-I0)+b4(ENt-E0)| .

(16)

這里,

SN(x,y,t),BN(x,y,t),IN(x,y,t)及EN(x,y,t)依次表示營養(yǎng)物濃度、活生物密度、死生物密度及胞外聚合物密度在點(diǎn)(x,y,t)處的數(shù)值解。

算例1 計(jì)算使用的參數(shù)μ=0.02,YH=0.003,YE=KE=0.03,KS=0.2,KI=0.04,KL=0.04;α=β=2,問題域?yàn)棣?(0,1),控制方程為(1)-(5),初始條件為:

S0(x)=1 ,

B0(x)=0.15e-200(x-0.3)2+0.2e-150(x-0.65)2,

I0(x)=0.04e-80(x-0.2)2+0.03e-60(x-0.7)2,

E0(x)=0.025 .

在邊界網(wǎng)格點(diǎn)上施加均勻的Neumann邊界條件，使用的網(wǎng)格為Δt=0.001 and Δx=0.01，數(shù)值結(jié)果如圖 1所示。

圖1 使用有限差分格式(6)-(9)得到的數(shù)值解

由式(16)計(jì)算出的5個(gè)時(shí)刻t=0.25,0.5,0.75,1,2的宏觀誤差依次為:

1.589 9×10-4,1.368 3×10-4,

1.131 9×10-4,8.583 0×10-5,3.162 4×10-4.

圖2給出的是4個(gè)不同時(shí)刻t=5,10,20,30的數(shù)值解。相應(yīng)的宏觀誤差依次為:2.685 6×10-5,2.180 2×10-6,9.097 0×10-5,1.178 6×10-4.在這種長期情況下,營養(yǎng)基質(zhì)函數(shù)是趨于零的,而活微生物擴(kuò)散到整個(gè)空間區(qū)域。對于長時(shí)間的行為,進(jìn)一步的實(shí)驗(yàn)證實(shí)活的微生物密度最終趨于零,這也是我們所期望的。

綜合圖 1和圖2,我們發(fā)現(xiàn),起初隨著營養(yǎng)基質(zhì)密度的減小,總的微生物量有所增加,隨著時(shí)間的增加,活性微生物轉(zhuǎn)化成惰性微生物。

從圖3可以觀察到,使用有限差分格式(11)-(14)得到的數(shù)值解具有奇異性和震蕩現(xiàn)象,說明該格式不穩(wěn)定,不能確保數(shù)值解是非負(fù)、有界的。而本文提出的有限差分格式是穩(wěn)定的,得到的數(shù)值解是非負(fù)、有界的。

算例2 計(jì)算使用的參數(shù)μ=2,YH=0.35,YE=KE=0.03,KS=0.2,KI=0.4,KL=0.001;α=β=4,問題域?yàn)棣?(0,1),控制方程為(1)-(5),初始條件為:

S0(x)=0.8，

B0(x)=0.15e-200(x-0.3)2+0.2e-150(x-0.65)2，

I0(x)=0.02e-80x2+0.01e-60(x-1)2，

E0(x)=0.005 .

在邊界網(wǎng)格點(diǎn)上施加均勻的Neumann邊界條件,使用的網(wǎng)格為Δt=0.001 and Δx=0.01,數(shù)值結(jié)果見圖4。

圖4給出了4個(gè)不同時(shí)刻t=1,2,3,10的數(shù)值解.相應(yīng)的宏觀誤差依次為:1.765 1×10-6,6.717 4×10-6,3.037 0×10-6,可以觀察到,在起初階段,隨著營養(yǎng)基質(zhì)函數(shù)的減小,活微生物密度總量增加;當(dāng)營養(yǎng)基質(zhì)趨于零時(shí),活微生物密度總減小,逐漸趨于零;而對惰性生物密度總量始終是增加的。這些結(jié)果符合系統(tǒng)的特征。

圖2 使用有限差分格式(6)-(9)得到的數(shù)值解

圖3 使用有限差分格式(11)-(14)得到的數(shù)值解。使用的控制方程、初始條件、邊界條件、計(jì)算參數(shù)以及計(jì)算網(wǎng)格和圖 1和圖2相同

圖4 使用有限差分格式(6)-(9)得到的數(shù)值解

從圖5可以觀察到,對數(shù)值算例2使用有限差分格式(11)-(14)得到的數(shù)值解也具有奇異性和震蕩現(xiàn)象,說明該格式不穩(wěn)定,不能確保數(shù)值解是非負(fù)、有界的.而本章提出的有限差分格式是穩(wěn)定的,得到的數(shù)值解是非負(fù)、有界的。

5 結(jié)論

本文用新的顯式有限差分法來逼近復(fù)雜生物膜系統(tǒng)的正的、有界解,其控制方程是由四個(gè)拋物線型偏微分方程構(gòu)成,偏微分方程中同樣也包含一個(gè)與密度相關(guān)的非線性擴(kuò)散反應(yīng)項(xiàng).理論上,模型的非負(fù)、有界解的存在性和唯一性文獻(xiàn)中已有證明,而獲得該模型的精確解確實(shí)是很困難的,因此,我們不得不采用數(shù)值模擬的方法來得到近似解。本文的理論和數(shù)值研究總結(jié)如下:

1) 建立了復(fù)雜生物膜系統(tǒng)的顯式有限差分格式;

2) 給出了2個(gè)數(shù)值算例，結(jié)果表明數(shù)值結(jié)果和實(shí)際是相吻合的，能夠反應(yīng)微生物成長規(guī)律;

3) 和標(biāo)準(zhǔn)中心差分方法相比較，本文所提出的方法具有較好的穩(wěn)定性，不會出現(xiàn)奇異的數(shù)值解和震蕩現(xiàn)象;

4) 不足之處，相關(guān)收斂性、穩(wěn)定性等理論體系還沒建立起來。

圖5 使用有限差分格式(11)-(14)得到的數(shù)值解。使用的控制方程、初始條件、邊界條件、計(jì)算參數(shù)以及計(jì)算網(wǎng)格和圖4相同

[1] 李治平.偏微分方程數(shù)值解講義[M].北京:北京大學(xué)出版社,2010.

[2] Liu WK,Liu YL,Farrell D.Immersed finite element method and its applications to biological systems[J].Computer Methods In Applied Mechanics And Engineering,2006,195(13-16):1722-1749.

[3] Liu G R,Quek S S.Finite Element Method:a practical course[M].Burlington MA:BH,2003.

[4] Wang Shuai,Yuan Guangwei,Li Yonghai.A monotone finite volume scheme for advection-diffusion equations on distorted meshes[J].International Journal for Numerical Methods in Fluids,2012,69(7):1283-1298.

[5] Hermeline F.A finite volume method for the approximation of convection-diffusion equations on general meshes[J].International Journal For Numerical Methods In Engineering,2012,91(12):1331-1357.

[6] Macías-Díaz J E.A bounded finite-difference discretization of a two-dimensional diffusion equation with logistic nonlinear reaction[J].Int J Mod Phys C,2011,22:953-966.

[7] Macías-Díaz J E,Jerez-Galiano S,Puri A.Positivity-preserving methods for a linearized Fisher-KPP equation with consistency properties in the energy domain[J].J Differ Equ Appl,2010,16:389-405.

[8] Macías-Díaz J E,Puri A.A boundedness-preserving finite-difference scheme for a dampednonlinear wave equation[J].Appl Numer Math,2010,60:934-948.

[9] Eberl H J,Parker D,van Loosdrecht M C M.A new deterministic spatio-temporal continuum model for biofilm development[J].J Theor Med,2001,3:161-176.

[10] Efendiev M A,Eberl H J,Zelik S V.Existence and longtime behavior ofsolutions of a nonlinear reaction-diffusion system arising in the modeling of biofilms[J].RIMS Kokyuroko,2002,1258:49-71.

[11] Eberl H J.A deterministic continuum model for the formation of EPS in heterogeneous boilarchitectures[J].Proc Biofilms,2004,237-242.

(編輯：朱 倩)

Explicit Finite-difference Method for a ComplexBio-film Dynamics System

SUN Gaofeng,NIU Ruiping,LI Ming

(CollegeofMathematics,TaiyuanUniversityofTechnology,Taiyuan030024,China)

Numerical methods with non-negative and bounded solutions have great significance because of the requirement from many mathematical models with practical background.In this paper,a non-negative and bounded explicit finite-difference scheme is presented for bio-film growth dynamics system.The scheme describes mutual process among different microorganism compositions on nutritional substance and is expressed by four nonlinear parabolic partial differential equations with variable coefficients.The explicit finite-difference scheme neatly transforms the nonlinear form into linear form,which ensures the non-negative and bounded property of solutions.Numerical examples are presented to illustrate the performance of the methods.Error in macro significance is also given.

complex biofilm dynamics system；explicit finite-difference method；nonlinear reaction-diffusion equations；Nonnegative bounded

1007-9432(2015)06-0790-07

2015-04-09

國家自然科學(xué)基金資助項(xiàng)目:可商業(yè)化光滑有限元分析軟件的研發(fā)(11472184)

孫高峰(1976-)， 男，山西臨猗人，博士，主要從事計(jì)算數(shù)學(xué)的研究，(Tel)13453157694

李明(1982-)，男，教授，主要從事無網(wǎng)絡(luò)計(jì)算方法，偏微方程中的反問題等的研究，(E-mail)liming01@tyut.edu.cn

O241

A

10.16355/j.cnki.issn1007-9432tyut.2015.06.028