反应扩散模型在图灵斑图中的应用及数值模拟∗

张荣培 王震 王语 韩子健

1)(沈阳师范大学数学与系统科学学院,沈阳 110034)

2)(山东科技大学数学与系统科学学院,青岛 266590)

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

1 引 言

斑图是在空间或时间上具有某种规律性的非均匀宏观结构,普遍存在于自然界.1952年,著名的英国数学家图灵把他的目光转向生物学领域,用一个反应扩散系统成功地说明了某些生物体表面图纹产生的原理[1].图灵从数学角度表明,在反应扩散系统中,稳定状态会在某些条件下失稳,并自发产生空间定态图纹,此斑图通常称为图灵斑图.

经过多年的研究,各界学者利用反应扩散系统预测得到了更多的图灵斑图,在理论和实验方面取得了许多重要成果.他们证实了化学系统中图灵斑图的形成[2],讨论自催化反应中的动力学行为,探讨此类耦合反应扩散体系中影响图灵斑图的因素[3].给出Gray-Scott模型、Brusselator模型等系统扩散引起不稳定的数学机理[4],并描述了Gierer-Meinhardt,Lengyel-Epstein等模型的某些动力学行为(性质)[5,6].最近几年,图灵斑图在实验方面取得一系列最新的进展,Copie等[7]运用实验在一个双稳态被动非线性共振器中探讨了图灵调制和法拉第参数不稳定性的相互作用;Tompkins等[8]利用微流体化学室证实图灵理论体系,并观测到第七种时空模式;Lacitignola[9]研究了图灵不稳定现象的发生条件,论述了具体形态的电化学反应扩散模型在一个球面上的图案形成的特性;Gaskins等[10]在二氧化氯碘丙二酸反应实验中,通过添加卤化钠盐溶液得到新的图灵斑图.

在这些系统中存在两种化学反应物质,它们不仅能相互作用,而且还能进行独自扩散.事实上,图灵斑图的产生对应的是一个非线性反应动力学过程与一种特殊扩散过程的耦合.这个特殊的扩散过程由于两种因子的扩散速度不同会发生失稳,这就是图灵斑图产生的机理.在数学上,图灵斑图可以用无量纲化的反应扩散方程组描述[11],即

式中u和v是系统变量,分别代表参与化学反应的两种物质的浓度;c和d是扩散系数,t是时间变量,f(u,v)和g(u,v)表示反应项.设Ω为RN中带有光滑边界的有界区域,Ω=[0,a]×[0,b],边界为∂Ω,边界条件为齐次Neumann边界条件,即其中n表示边界上单位外法向.

由于(1)式为耦合的非线性反应扩散方程,很难得到其精确解.近年来,许多学者用有限差分方法、有限元方法、谱方法等[12−14]多种数值方法求解(1)式,这些方法各有特点.相比于有限元方法和有限差分方法的低阶精度,谱方法[14]仅用少量的节点,采用Legendre,Chebyshev等适合的正交多项式离散即可达到指数阶收敛的谱精度.图灵斑图在空间上的结构具有一定的规律,且解比较光滑,因此采用谱方法离散是可行的.常用的谱配置方法主要有Fourier配置法[15],Chebyshev配置法[16],Hermite配置法等[17].由于本文考虑的(1)式边界条件为齐次Neumann边界条件,因此采用Chebyshev配置方法求解(1)式.

对(1)式进行空间离散后,得到的是刚性的非线性常微分方程组(ODEs).显式时间离散方法虽可以用迭代的方法求解,但其对时间步长有严格的约束;隐式方法虽然可以允许大的时间步长,但是对于阶数非常大的非线性方程组的求解问题十分复杂,这对于全隐式方法来说是一个巨大的挑战.由于谱配置法所得到的谱微分矩阵是满的,显然利用追赶法等代数线性方程组的快速解法是不合适的,因此交替方向隐式方法在这里并不适用.本文采用紧致隐积分因子(compact implicit integration factor,cIIF)方法求解ODEs.2006年Nie等[18]以隐积分因子(IIF)方法为基础发展了cIIF方法.传统的隐积分因子方法在求解高维问题时,离散矩阵的指数运算的存储量和运算量非常大,导致运算速度缓慢.紧致隐积分因子方法[19]通过引入离散矩阵的紧致表达式并在各个方向进行矩阵的指数运算,使得中央处理器(CPU)的存储大大降低,计算速度也得到了显著提高.

本文内容安排如下:第2节对反应扩散方程组进行线性分析,通过特征值解释图灵斑图的数学机理,然后以Gierer-Meinhardt模型为例分析系统处于稳定状态和不稳定状态时各参数需要满足的条件,进而探索斑图形成需要满足的条件;第3节研究数值方法,在空间离散条件下采用Chebyshev谱方法,时间离散条件下采用紧致隐积分因子方法,用MATLAB进行编程求解;第4节给出大量数值实验并对理论分析结果进行验证.

2 图灵斑图的形成

2.1 斑图形成的数学机理

首先考虑(1)式没有扩散项,假设存在惟一的均匀定态解(u0,v0),即常数u0,v0满足

令U=u−u0,V=v−v0,并在(u0,v0)处线性化后得到如下系统:

式中c11=fu(u0,v0),c12=fv(u0,v0),c21=gu(u0,v0),c22=gv(u0,v0).均匀定态解(u0,v0)在没有扩散时是稳定的,这等价于相应的特征值问题的矩阵的特征值实部是负数.

考虑加入扩散项后的反应扩散方程组((1)式).如果此时产生斑图,即(u0,v0)是不稳定的,要求特征值有正实部.所谓不稳定,体现为两种反应物的扩散速度不同,从而引起失稳.对(1)式作线性化处理,研究特征值正实部引起的线性不稳定性,进而推导出原方程的不稳定性.对均匀定态解(u0,v0)作一个微扰,可得线性微扰方程为

求解如下方程可得相应的特征值:

式中λ为特征值.只要(5)式中的特征值有正实部,则(u0,v0)对于(1)式是不稳定的.考虑到齐次Neumann边界条件,得到(5)式所对应的特征值为

具体推导过程见附录A.

2.2 Gierer-Meinhardt模型

生物的发育过程是复杂的,其中重要的是形态形成阶段,与之对应的是生物体内器官的形成.由于该阶段的重要性,渐渐形成一个新的领域——形态学,主要研究导致细胞分化和定位因素的浓度对组织器官的影响.Gierer-Meinhardt模型是由Gierer和Meinhardt在研究激活物和抑制剂两种不同物质的产生和扩散时建立的[20],之后Gierer和Meinhardt利用数值方法导出一维和二维空间区域中上述系统产生多样斑图的条件.Gierer-Meinhardt模型被广泛应用于形态形成过程中一些基本现象的研究,最近的一些工作可以参见文献[21—23].

以Gierer-Meinhardt模型为例,结合上述理论分析,计算产生斑图时需要满足的条件.取(1)式中

其中系数κ,η,ε为系统的控制参数,固定η=0.1,ε=0.04.由此得到线性化系统(3)式中的系数为

易得该系统的特征值为λ1=−1.2984,λ2=−7.7016,此时系统是稳定的.

加入扩散项后,原方程组对应的特征问题为

相应的特征方程为

为使(8)式含有正实部的特征值,需要考虑两种情况.

第一种情况是两个特征值异号,则应满足

图1 特征值的实部Re(λ)随参数的变化 (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.

第二种情况是两个特征值都是正的,应满足此时κ无解.

由于反应扩散方程组联系于解析半群,所以线性化后的正实部特征值引起的不稳定性可以推导出原方程组的不稳定性.故当κ＞κ0=0.0093248时,系统处于不稳定状态,因而系统能够产生斑图.特征值的实部Re(λ)在参数κ取不同值时的变化如图1所示.

从图1可以看出,当κ= 0.0128＞κ0和κ=0.0152＞κ0时,特征值的实部会出现正值,此时系统不稳定;当κ=0.008＜κ0时,特征值的实部始终为负,系统最后会达到稳定状态.第3节将用数值算例验证该结论.

3 数值方法

3.1 Chebyshev谱配置法

将求解区域[−1,1]2离散为Gauss-Lobatto网格,即

其中Nx和Ny是正整数.对于一般的求解区域Ω=[a,b]×[c,d],可以采用公式

将区域转化为[−1,1]2.在网格Th中将u(x,y)数值解定义为矩阵形式,U∈R(Nx−1)×(Ny−1),

式中ui,j表示u在网格点(xi,xj)的数值解.引入Chebyshev一阶微分矩阵和二阶微分矩阵(具体推导过程见附录B).则u(x,y)关于x的二阶偏导数在配置点的值,可以用矩阵乘积的形式近似为矩阵Ax是在Chebyshev二阶微分矩阵基础上考虑Neumann边界条件得到的,

其中

同样地,对于y的二阶偏导数,有UAy,其中矩阵Ay定义同Ax.借助谱微分矩阵,可将方程中的Laplace算子离散成矩阵乘积的形式,即将Chebyshev谱配置方法应用于反应扩散方程,得到其半离散形式为

3.2 紧致隐积分因子法

将对空间离散后得到的非线性常微分方程组((12)式)采用紧致隐积分因子方法进行时间离散.定义时间步长为τ=Δt,第n层时间步为tn=nτ,n=0,1,2,···. 在(12)式两端同时左乘指数矩阵e−Axt,右乘指数矩阵e−Ayt.为描述方便,取(10)式c=1,d=1,可将(12)式中第一个等式写为

将时间离散为0=t0＜t1＜···,将(13)式在一个时间步长内关于时间积分,并用梯形公式近似可得二阶紧致隐积分因子格式为

进一步化简得

在非线性方程组(13)式中,右端第一项可以通过矩阵乘积得到,右端第二项采用Picard迭代方法求解:

同理处理(12)式中第二个等式可得

该方法中矩阵eAxΔt和eAyΔt的阶数分别为Nx×Nx和Ny×Ny.在空间网格剖分量很大时,该方法可以降低存储量和运算量,使计算速度更快.

4 数值算例

对于前述Gierer-Meinhardt模型,取Ω=(−1,1)×(−1,1),η=0.1,c=0.04,κ是不固定的参数.设

其中

图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 数值算例I

取κ=0.0128,N=100,h=2/100=0.02,τ=0.1h,t取图2所示各值时,得到对应的图像.由图2可知,随着时间的推移,初始扰动不断增强扩大,最终形成清晰的斑图.

4.2 数值算例II

取κ=0.0152,t取图3所示各值,其他参数与算例I相同,可得到t取不同值时对应的图像.由图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 数值算例III

取κ=0.008,其他取值与算例II相同,t取不同值时对应的图像如图4所示.由图4可知,随着时间的推移,系统达到稳定状态,反应扩散模型不能形成斑图.

由数值模拟结果来看,其他条件一定的情况下,κ取不同值对于产生斑图有重要的影响.数值模拟结果与理论结果一致.

此外,我们也对周期性边界条件的Gierer-Meinhardt模型采用Fourier谱方法进行数值求解,结果显示周期边界条件对斑图的形状几乎没有影响.

5 结 论

介绍了图灵斑图形成的数学机理,并结合Gierer-Meinhardt模型,分析系统不稳定状态的各系数需要满足的条件,即产生斑图的条件.运用紧致隐积分因子方法大大减少了存储和CPU运算时间,该方法对于大时间数值模拟是一个高效、高精度的数值方法.数值算例模拟了斑图形成的过程,验证了理论分析结果.这些结论还可应用于求解带有分数阶的反应扩散方程组.

图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 图灵斑图的形成机理

首先在区域Ω⊂RN(N=1,2)内考虑带有齐次Neumann边界条件的Laplace算子的特征值问题.一维情况下,特征值问题为

式中a∈R+.特征值问题可表示为µ2−λ=0,解得只有λ＜0时可解得特征值λk=−(kπ/a)2,且特征值所对应的特征函数为

在二维情况下,特征值问题为

式中a,b∈R+,应采用分离变量法求解特征值. 设u=X(x)Y(y),代入方程得设解得故特征值为λk,l=且特征值所对应的特征函数为

考虑方程组的特征值问题,令

代入原方程组可得

当方程组(A3)有非零解,满足

此时方程组所对应的特征值为

附录B 谱微分矩阵

定义在[−1,1]上的标准k阶Chebyshev多项式Tk(x)为Tk(x)=cos(karccosx),k=0,1,2,···. 令x=cosz,则有Tk=coskz,满足如下递推关系:

Tk(x)在[−1,1]上的N+1个Gauss-Lobatto点值为零:

设N阶多项式uN(x)∈PN在上述配置点xj满足uN(xj)=u(xj),则有

式中hj(x)为N阶Lagrange基函数.用配置法求解未知量在网格点处的值,需要表示配置点处的导数值.对(B3)式求p阶导数,得

[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物理学报62 220503]

[3]Zhang L,Liu S Y 2007Appl.Math.Mec.28 1102(in Chinese)[张丽,刘三阳2007应用数学和力学28 1102]

[4]Li B,Wang M X 2008Appl.Math.Mec.29 749(in Chinese)[李波,王明新2008应用数学和力学29 749]

[5]Hu W Y,Shao Y Z 2014Acta Phys.Sin.63 238202(in Chinese)[胡文勇,邵元智 2014物理学报 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物理学报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东北师大学报3 26]