張榮培 王震 王語 韓子健

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

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

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

1 引 言

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

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

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

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

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

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

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

2 圖靈斑圖的形成

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

首先考慮(1)式?jī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)在沒有擴散時是穩(wěn)定的,這等價于相應(yīng)的特征值問題的矩陣的特征值實部是負數(shù).

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

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

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

具體推導(dǎo)過程見附錄A.

2.2 Gierer-Meinhardt模型

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

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

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

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

加入擴散項后,原方程組對應(yīng)的特征問題為

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

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

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

圖1 特征值的實部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.

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

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

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

3 數(shù)值方法

3.1 Chebyshev譜配置法

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

其中Nx和Ny是正整數(shù).對于一般的求解區(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)格點(xi,xj)的數(shù)值解.引入Chebyshev一階微分矩陣和二階微分矩陣(具體推導(dǎo)過程見附錄B).則u(x,y)關(guān)于x的二階偏導(dǎo)數(shù)在配置點的值,可以用矩陣乘積的形式近似為矩陣Ax是在Chebyshev二階微分矩陣基礎(chǔ)上考慮Neumann邊界條件得到的,

其中

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

3.2 緊致隱積分因子法

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

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

進一步化簡得

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

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

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

4 數(shù)值算例

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

其中

圖2 取κ=0.0128時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所示各值時,得到對應(yīng)的圖像.由圖2可知,隨著時間的推移,初始擾動不斷增強擴大,最終形成清晰的斑圖.

4.2 數(shù)值算例II

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

圖3 取κ=0.0152時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取不同值時對應(yīng)的圖像如圖4所示.由圖4可知,隨著時間的推移,系統(tǒng)達到穩(wěn)定狀態(tài),反應(yīng)擴散模型不能形成斑圖.

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

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

5 結(jié) 論

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

圖4 取κ=0.008時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 圖靈斑圖的形成機理

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

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

在二維情況下,特征值問題為

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

考慮方程組的特征值問題,令

代入原方程組可得

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

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

附錄B 譜微分矩陣

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

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

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

式中hj(x)為N階Lagrange基函數(shù).用配置法求解未知量在網(wǎng)格點處的值,需要表示配置點處的導(dǎo)數(shù)值.對(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)[李新政,白占國,李燕,趙昆,賀亞峰2013物理學(xué)報62 220503]

[3]Zhang L,Liu S Y 2007Appl.Math.Mec.28 1102(in Chinese)[張麗,劉三陽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é)報 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)[白占國,董麗芳,李永輝,范偉麗2011物理學(xué)報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é)報3 26]