反應(yīng)擴(kuò)散模型在圖靈斑圖中的應(yīng)用及數(shù)值模擬?

張榮培 王震 王語(yǔ) 韓子健

1)(沈陽(yáng)師范大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,沈陽(yáng) 110034)

2)(山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,青島 266590)

(2017年8月6日收到;2017年11月6日收到修改稿)

1 引 言

斑圖是在空間或時(shí)間上具有某種規(guī)律性的非均勻宏觀結(jié)構(gòu),普遍存在于自然界.1952年,著名的英國(guó)數(shù)學(xué)家圖靈把他的目光轉(zhuǎn)向生物學(xué)領(lǐng)域,用一個(gè)反應(yīng)擴(kuò)散系統(tǒng)成功地說(shuō)明了某些生物體表面圖紋產(chǎn)生的原理[1].圖靈從數(shù)學(xué)角度表明,在反應(yīng)擴(kuò)散系統(tǒng)中,穩(wěn)定狀態(tài)會(huì)在某些條件下失穩(wěn),并自發(fā)產(chǎn)生空間定態(tài)圖紋,此斑圖通常稱為圖靈斑圖.

經(jīng)過(guò)多年的研究,各界學(xué)者利用反應(yīng)擴(kuò)散系統(tǒng)預(yù)測(cè)得到了更多的圖靈斑圖,在理論和實(shí)驗(yàn)方面取得了許多重要成果.他們證實(shí)了化學(xué)系統(tǒng)中圖靈斑圖的形成[2],討論自催化反應(yīng)中的動(dòng)力學(xué)行為,探討此類耦合反應(yīng)擴(kuò)散體系中影響圖靈斑圖的因素[3].給出Gray-Scott模型、Brusselator模型等系統(tǒng)擴(kuò)散引起不穩(wěn)定的數(shù)學(xué)機(jī)理[4],并描述了Gierer-Meinhardt,Lengyel-Epstein等模型的某些動(dòng)力學(xué)行為(性質(zhì))[5,6].最近幾年,圖靈斑圖在實(shí)驗(yàn)方面取得一系列最新的進(jìn)展,Copie等[7]運(yùn)用實(shí)驗(yàn)在一個(gè)雙穩(wěn)態(tài)被動(dòng)非線性共振器中探討了圖靈調(diào)制和法拉第參數(shù)不穩(wěn)定性的相互作用;Tompkins等[8]利用微流體化學(xué)室證實(shí)圖靈理論體系,并觀測(cè)到第七種時(shí)空模式;Lacitignola[9]研究了圖靈不穩(wěn)定現(xiàn)象的發(fā)生條件,論述了具體形態(tài)的電化學(xué)反應(yīng)擴(kuò)散模型在一個(gè)球面上的圖案形成的特性;Gaskins等[10]在二氧化氯碘丙二酸反應(yīng)實(shí)驗(yàn)中,通過(guò)添加鹵化鈉鹽溶液得到新的圖靈斑圖.

在這些系統(tǒng)中存在兩種化學(xué)反應(yīng)物質(zhì),它們不僅能相互作用,而且還能進(jìn)行獨(dú)自擴(kuò)散.事實(shí)上,圖靈斑圖的產(chǎn)生對(duì)應(yīng)的是一個(gè)非線性反應(yīng)動(dòng)力學(xué)過(guò)程與一種特殊擴(kuò)散過(guò)程的耦合.這個(gè)特殊的擴(kuò)散過(guò)程由于兩種因子的擴(kuò)散速度不同會(huì)發(fā)生失穩(wěn),這就是圖靈斑圖產(chǎn)生的機(jī)理.在數(shù)學(xué)上,圖靈斑圖可以用無(wú)量綱化的反應(yīng)擴(kuò)散方程組描述[11],即

式中u和v是系統(tǒng)變量,分別代表參與化學(xué)反應(yīng)的兩種物質(zhì)的濃度;c和d是擴(kuò)散系數(shù),t是時(shí)間變量,f(u,v)和g(u,v)表示反應(yīng)項(xiàng).設(shè)?為RN中帶有光滑邊界的有界區(qū)域,?=[0,a]×[0,b],邊界為??,邊界條件為齊次Neumann邊界條件,即其中n表示邊界上單位外法向.

由于(1)式為耦合的非線性反應(yīng)擴(kuò)散方程,很難得到其精確解.近年來(lái),許多學(xué)者用有限差分方法、有限元方法、譜方法等[12?14]多種數(shù)值方法求解(1)式,這些方法各有特點(diǎn).相比于有限元方法和有限差分方法的低階精度,譜方法[14]僅用少量的節(jié)點(diǎn),采用Legendre,Chebyshev等適合的正交多項(xiàng)式離散即可達(dá)到指數(shù)階收斂的譜精度.圖靈斑圖在空間上的結(jié)構(gòu)具有一定的規(guī)律,且解比較光滑,因此采用譜方法離散是可行的.常用的譜配置方法主要有Fourier配置法[15],Chebyshev配置法[16],Hermite配置法等[17].由于本文考慮的(1)式邊界條件為齊次Neumann邊界條件,因此采用Chebyshev配置方法求解(1)式.

對(duì)(1)式進(jìn)行空間離散后,得到的是剛性的非線性常微分方程組(ODEs).顯式時(shí)間離散方法雖可以用迭代的方法求解,但其對(duì)時(shí)間步長(zhǎng)有嚴(yán)格的約束;隱式方法雖然可以允許大的時(shí)間步長(zhǎng),但是對(duì)于階數(shù)非常大的非線性方程組的求解問(wèn)題十分復(fù)雜,這對(duì)于全隱式方法來(lái)說(shuō)是一個(gè)巨大的挑戰(zhàn).由于譜配置法所得到的譜微分矩陣是滿的,顯然利用追趕法等代數(shù)線性方程組的快速解法是不合適的,因此交替方向隱式方法在這里并不適用.本文采用緊致隱積分因子(compact implicit integration factor,cIIF)方法求解ODEs.2006年Nie等[18]以隱積分因子(IIF)方法為基礎(chǔ)發(fā)展了cIIF方法.傳統(tǒng)的隱積分因子方法在求解高維問(wèn)題時(shí),離散矩陣的指數(shù)運(yùn)算的存儲(chǔ)量和運(yùn)算量非常大,導(dǎo)致運(yùn)算速度緩慢.緊致隱積分因子方法[19]通過(guò)引入離散矩陣的緊致表達(dá)式并在各個(gè)方向進(jìn)行矩陣的指數(shù)運(yùn)算,使得中央處理器(CPU)的存儲(chǔ)大大降低,計(jì)算速度也得到了顯著提高.

本文內(nèi)容安排如下:第2節(jié)對(duì)反應(yīng)擴(kuò)散方程組進(jìn)行線性分析,通過(guò)特征值解釋圖靈斑圖的數(shù)學(xué)機(jī)理,然后以Gierer-Meinhardt模型為例分析系統(tǒng)處于穩(wěn)定狀態(tài)和不穩(wěn)定狀態(tài)時(shí)各參數(shù)需要滿足的條件,進(jìn)而探索斑圖形成需要滿足的條件;第3節(jié)研究數(shù)值方法,在空間離散條件下采用Chebyshev譜方法,時(shí)間離散條件下采用緊致隱積分因子方法,用MATLAB進(jìn)行編程求解;第4節(jié)給出大量數(shù)值實(shí)驗(yàn)并對(duì)理論分析結(jié)果進(jìn)行驗(yàn)證.

2 圖靈斑圖的形成

2.1 斑圖形成的數(shù)學(xué)機(jī)理

首先考慮(1)式?jīng)]有擴(kuò)散項(xiàng),假設(shè)存在惟一的均勻定態(tài)解(u0,v0),即常數(shù)u0,v0滿足

令U=u?u0,V=v?v0,并在(u0,v0)處線性化后得到如下系統(tǒng):

式中c11=fu(u0,v0),c12=fv(u0,v0),c21=gu(u0,v0),c22=gv(u0,v0).均勻定態(tài)解(u0,v0)在沒(méi)有擴(kuò)散時(shí)是穩(wěn)定的,這等價(jià)于相應(yīng)的特征值問(wèn)題的矩陣的特征值實(shí)部是負(fù)數(shù).

考慮加入擴(kuò)散項(xiàng)后的反應(yīng)擴(kuò)散方程組((1)式).如果此時(shí)產(chǎn)生斑圖,即(u0,v0)是不穩(wěn)定的,要求特征值有正實(shí)部.所謂不穩(wěn)定,體現(xiàn)為兩種反應(yīng)物的擴(kuò)散速度不同,從而引起失穩(wěn).對(duì)(1)式作線性化處理,研究特征值正實(shí)部引起的線性不穩(wěn)定性,進(jìn)而推導(dǎo)出原方程的不穩(wěn)定性.對(duì)均勻定態(tài)解(u0,v0)作一個(gè)微擾,可得線性微擾方程為

求解如下方程可得相應(yīng)的特征值:

式中λ為特征值.只要(5)式中的特征值有正實(shí)部,則(u0,v0)對(duì)于(1)式是不穩(wěn)定的.考慮到齊次Neumann邊界條件,得到(5)式所對(duì)應(yīng)的特征值為

具體推導(dǎo)過(guò)程見(jiàn)附錄A.

2.2 Gierer-Meinhardt模型

生物的發(fā)育過(guò)程是復(fù)雜的,其中重要的是形態(tài)形成階段,與之對(duì)應(yīng)的是生物體內(nèi)器官的形成.由于該階段的重要性,漸漸形成一個(gè)新的領(lǐng)域——形態(tài)學(xué),主要研究導(dǎo)致細(xì)胞分化和定位因素的濃度對(duì)組織器官的影響.Gierer-Meinhardt模型是由Gierer和Meinhardt在研究激活物和抑制劑兩種不同物質(zhì)的產(chǎn)生和擴(kuò)散時(shí)建立的[20],之后Gierer和Meinhardt利用數(shù)值方法導(dǎo)出一維和二維空間區(qū)域中上述系統(tǒng)產(chǎn)生多樣斑圖的條件.Gierer-Meinhardt模型被廣泛應(yīng)用于形態(tài)形成過(guò)程中一些基本現(xiàn)象的研究,最近的一些工作可以參見(jiàn)文獻(xiàn)[21—23].

以Gierer-Meinhardt模型為例,結(jié)合上述理論分析,計(jì)算產(chǎn)生斑圖時(shí)需要滿足的條件.取(1)式中

其中系數(shù)κ,η,ε為系統(tǒng)的控制參數(shù),固定η=0.1,ε=0.04.由此得到線性化系統(tǒng)(3)式中的系數(shù)為

易得該系統(tǒng)的特征值為λ1=?1.2984,λ2=?7.7016,此時(shí)系統(tǒng)是穩(wěn)定的.

加入擴(kuò)散項(xiàng)后,原方程組對(duì)應(yīng)的特征問(wèn)題為

相應(yīng)的特征方程為

為使(8)式含有正實(shí)部的特征值,需要考慮兩種情況.

第一種情況是兩個(gè)特征值異號(hào),則應(yīng)滿足

圖1 特征值的實(shí)部Re(λ)隨參數(shù)的變化 (a)κ=0.0128;(b)κ=0.0152;(c)κ=0.008Fig.1.Real part Re(λ)of eigenvalues varying with parameters:(a)κ=0.0128;(b)κ=0.0152;(c)κ=0.008.

第二種情況是兩個(gè)特征值都是正的,應(yīng)滿足此時(shí)κ無(wú)解.

由于反應(yīng)擴(kuò)散方程組聯(lián)系于解析半群,所以線性化后的正實(shí)部特征值引起的不穩(wěn)定性可以推導(dǎo)出原方程組的不穩(wěn)定性.故當(dāng)κ＞κ0=0.0093248時(shí),系統(tǒng)處于不穩(wěn)定狀態(tài),因而系統(tǒng)能夠產(chǎn)生斑圖.特征值的實(shí)部Re(λ)在參數(shù)κ取不同值時(shí)的變化如圖1所示.

從圖1可以看出,當(dāng)κ= 0.0128＞κ0和κ=0.0152＞κ0時(shí),特征值的實(shí)部會(huì)出現(xiàn)正值,此時(shí)系統(tǒng)不穩(wěn)定;當(dāng)κ=0.008＜κ0時(shí),特征值的實(shí)部始終為負(fù),系統(tǒng)最后會(huì)達(dá)到穩(wěn)定狀態(tài).第3節(jié)將用數(shù)值算例驗(yàn)證該結(jié)論.

3 數(shù)值方法

3.1 Chebyshev譜配置法

將求解區(qū)域[?1,1]2離散為Gauss-Lobatto網(wǎng)格,即

其中Nx和Ny是正整數(shù).對(duì)于一般的求解區(qū)域?=[a,b]×[c,d],可以采用公式

將區(qū)域轉(zhuǎn)化為[?1,1]2.在網(wǎng)格Th中將u(x,y)數(shù)值解定義為矩陣形式,U∈R(Nx?1)×(Ny?1),

式中ui,j表示u在網(wǎng)格點(diǎn)(xi,xj)的數(shù)值解.引入Chebyshev一階微分矩陣和二階微分矩陣(具體推導(dǎo)過(guò)程見(jiàn)附錄B).則u(x,y)關(guān)于x的二階偏導(dǎo)數(shù)在配置點(diǎn)的值,可以用矩陣乘積的形式近似為矩陣Ax是在Chebyshev二階微分矩陣基礎(chǔ)上考慮Neumann邊界條件得到的,

其中

同樣地,對(duì)于y的二階偏導(dǎo)數(shù),有UAy,其中矩陣Ay定義同Ax.借助譜微分矩陣,可將方程中的Laplace算子離散成矩陣乘積的形式,即將Chebyshev譜配置方法應(yīng)用于反應(yīng)擴(kuò)散方程,得到其半離散形式為

3.2 緊致隱積分因子法

將對(duì)空間離散后得到的非線性常微分方程組((12)式)采用緊致隱積分因子方法進(jìn)行時(shí)間離散.定義時(shí)間步長(zhǎng)為τ=Δt,第n層時(shí)間步為tn=nτ,n=0,1,2,···. 在(12)式兩端同時(shí)左乘指數(shù)矩陣e?Axt,右乘指數(shù)矩陣e?Ayt.為描述方便,取(10)式c=1,d=1,可將(12)式中第一個(gè)等式寫為

將時(shí)間離散為0=t0＜t1＜···,將(13)式在一個(gè)時(shí)間步長(zhǎng)內(nèi)關(guān)于時(shí)間積分,并用梯形公式近似可得二階緊致隱積分因子格式為

進(jìn)一步化簡(jiǎn)得

在非線性方程組(13)式中,右端第一項(xiàng)可以通過(guò)矩陣乘積得到,右端第二項(xiàng)采用Picard迭代方法求解:

同理處理(12)式中第二個(gè)等式可得

該方法中矩陣eAxΔt和eAyΔt的階數(shù)分別為Nx×Nx和Ny×Ny.在空間網(wǎng)格剖分量很大時(shí),該方法可以降低存儲(chǔ)量和運(yùn)算量,使計(jì)算速度更快.

4 數(shù)值算例

對(duì)于前述Gierer-Meinhardt模型,取?=(?1,1)×(?1,1),η=0.1,c=0.04,κ是不固定的參數(shù).設(shè)

其中

圖2 取κ=0.0128時(shí)Gierer-Meinhardt模型形成的斑圖 (a)t=20;(b)t=80;(c)t=170;(d)t=270;(e)t=320;(f)t=340;(g)t=500;(h)t=600;(i)t=900Fig.2.Turing patterns in Gierer-Meinhardt model when κ=0.0128:(a)t=20;(b)t=80;(c)t=170;(d)t=270;(e)t=320;(f)t=340;(g)t=500;(h)t=600;(i)t=900.

4.1 數(shù)值算例I

取κ=0.0128,N=100,h=2/100=0.02,τ=0.1h,t取圖2所示各值時(shí),得到對(duì)應(yīng)的圖像.由圖2可知,隨著時(shí)間的推移,初始擾動(dòng)不斷增強(qiáng)擴(kuò)大,最終形成清晰的斑圖.

4.2 數(shù)值算例II

取κ=0.0152,t取圖3所示各值,其他參數(shù)與算例I相同,可得到t取不同值時(shí)對(duì)應(yīng)的圖像.由圖3可知,隨著時(shí)間的推移,初始擾動(dòng)不斷增強(qiáng)擴(kuò)大,最終形成清晰的斑圖.

圖3 取κ=0.0152時(shí)Gierer-Meinhardt模型形成的斑圖 (a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=290;(g)t=520;(h)t=620;(i)t=990Fig.3.Turing patterns in Gierer-Meinhardt model when κ=0.0152:(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=290;(g)t=520;(h)t=620;(i)t=990.

4.3 數(shù)值算例III

取κ=0.008,其他取值與算例II相同,t取不同值時(shí)對(duì)應(yīng)的圖像如圖4所示.由圖4可知,隨著時(shí)間的推移,系統(tǒng)達(dá)到穩(wěn)定狀態(tài),反應(yīng)擴(kuò)散模型不能形成斑圖.

由數(shù)值模擬結(jié)果來(lái)看,其他條件一定的情況下,κ取不同值對(duì)于產(chǎn)生斑圖有重要的影響.數(shù)值模擬結(jié)果與理論結(jié)果一致.

此外,我們也對(duì)周期性邊界條件的Gierer-Meinhardt模型采用Fourier譜方法進(jìn)行數(shù)值求解,結(jié)果顯示周期邊界條件對(duì)斑圖的形狀幾乎沒(méi)有影響.

5 結(jié) 論

介紹了圖靈斑圖形成的數(shù)學(xué)機(jī)理,并結(jié)合Gierer-Meinhardt模型,分析系統(tǒng)不穩(wěn)定狀態(tài)的各系數(shù)需要滿足的條件,即產(chǎn)生斑圖的條件.運(yùn)用緊致隱積分因子方法大大減少了存儲(chǔ)和CPU運(yùn)算時(shí)間,該方法對(duì)于大時(shí)間數(shù)值模擬是一個(gè)高效、高精度的數(shù)值方法.數(shù)值算例模擬了斑圖形成的過(guò)程,驗(yàn)證了理論分析結(jié)果.這些結(jié)論還可應(yīng)用于求解帶有分?jǐn)?shù)階的反應(yīng)擴(kuò)散方程組.

圖4 取κ=0.008時(shí)Gierer-Meinhardt模型形成的斑圖 (a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=220;(g)t=290;(h)t=270;(i)t=990Fig.4.Turing patterns in Gierer-Meinhardt model when κ=0.008:(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=220;(g)t=290;(h)t=270;(i)t=990.

附錄A 圖靈斑圖的形成機(jī)理

首先在區(qū)域??RN(N=1,2)內(nèi)考慮帶有齊次Neumann邊界條件的Laplace算子的特征值問(wèn)題.一維情況下,特征值問(wèn)題為

式中a∈R+.特征值問(wèn)題可表示為μ2?λ=0,解得只有λ＜0時(shí)可解得特征值λk=?(kπ/a)2,且特征值所對(duì)應(yīng)的特征函數(shù)為

在二維情況下,特征值問(wèn)題為

式中a,b∈R+,應(yīng)采用分離變量法求解特征值. 設(shè)u=X(x)Y(y),代入方程得設(shè)解得故特征值為λk,l=且特征值所對(duì)應(yīng)的特征函數(shù)為

考慮方程組的特征值問(wèn)題,令

代入原方程組可得

當(dāng)方程組(A3)有非零解,滿足

此時(shí)方程組所對(duì)應(yīng)的特征值為

附錄B 譜微分矩陣

定義在[?1,1]上的標(biāo)準(zhǔn)k階Chebyshev多項(xiàng)式Tk(x)為Tk(x)=cos(karccosx),k=0,1,2,···. 令x=cosz,則有Tk=coskz,滿足如下遞推關(guān)系:

Tk(x)在[?1,1]上的N+1個(gè)Gauss-Lobatto點(diǎn)值為零:

設(shè)N階多項(xiàng)式uN(x)∈PN在上述配置點(diǎn)xj滿足uN(xj)=u(xj),則有

式中hj(x)為N階Lagrange基函數(shù).用配置法求解未知量在網(wǎng)格點(diǎn)處的值,需要表示配置點(diǎn)處的導(dǎo)數(shù)值.對(duì)(B3)式求p階導(dǎo)數(shù),得

[1]Turing A M 1952Philos.Trans.R.Soc.Lond.B2 37

[2]Li X Z,Bai Z G,Li Y,Zhao K,He Y F 2013Acta Phys.Sin.62 220503(in Chinese)[李新政,白占國(guó),李燕,趙昆,賀亞峰2013物理學(xué)報(bào)62 220503]

[3]Zhang L,Liu S Y 2007Appl.Math.Mec.28 1102(in Chinese)[張麗,劉三陽(yáng)2007應(yīng)用數(shù)學(xué)和力學(xué)28 1102]

[4]Li B,Wang M X 2008Appl.Math.Mec.29 749(in Chinese)[李波,王明新2008應(yīng)用數(shù)學(xué)和力學(xué)29 749]

[5]Hu W Y,Shao Y Z 2014Acta Phys.Sin.63 238202(in Chinese)[胡文勇,邵元智 2014物理學(xué)報(bào) 63 238202]

[6]Peng R Wang M 2007Sci.China A50 377

[7]Copie F,Conforti M,Kudlinski A,Mussot A,Trillo S 2016Phys.Rev.Lett.116 143901

[8]Tompkins N,Li N,Girabawe C,Heymann M,Ermentrout G B,Epstein I R,Fraden S 2014Proc.Natl.Acad.Sci.USA111 4397

[9]Lacitignola D,Bozzini B,Frittelli M,Sgura I 2017Commun.Nonlinear Sci.Numer.Simul.48 484

[10]Gaskins D K,Pruc E E,Epstein I R,Dolnik M 2016Phys.Rev.Lett.117 056001

[11]Zhang R P,Yu X J,Zhu J,Loula A 2014Appl.Math.Model.38 1612

[12]Zhang R P,Zhu J,Loula A,Yu X J 2016J.Comput.Appl.Math.302 312

[13]Bai Z G,Dong L F,Li Y H,Fan W L 2011Acta Phys.Sin.60 118201(in Chinese)[白占國(guó),董麗芳,李永輝,范偉麗2011物理學(xué)報(bào)60 118201]

[14]Zhang R,Zhu J,Yu X,Li M,Loula A F D 2017Appl.Math.Comput.310 194

[15]Lv Z Q,Zhang L M,Wang Y S 2014Chin.Phys.B23 120203

[16]Wang H 2010Comput.Phys.Commun.181 325

[17]Hoz F D L,Vadillo F 2013Commun.Comput.Phys.14 1001

[18]Nie Q,Zhang Y T,Zhao R 2006J.Comput.Phys.214 521

[19]Nie Q,Wan F Y M,Zhang Y T,Liu X F 2008J.Comput.Phys.227 5238

[20]Gierer A,Meinhardt H 1972Kybernetik12 30

[21]Ward M J,Wei J 2003J.Nonlinear Sci.13 209

[22]Wei J,Winter M 2004J.Math.Pures Appl.83 433

[23]Li H X 2015J.Northeast Normal University3 26(in Chinese)[李海俠2015東北師大學(xué)報(bào)3 26]