斑图动力学中Duffet-Boissonade 方程的数值模拟

彭芳麟，鞠筱林，刘 畅

(北京师范大学 物理学系，北京 100875)



斑图动力学中Duffet-Boissonade 方程的数值模拟

彭芳麟，鞠筱林，刘 畅

(北京师范大学 物理学系，北京 100875)

对斑图动力学中的Duffet-Boissonade反应扩散方程组作了一维和二维的数值模拟.设置了4种不同的初始条件并按照图灵分岔的条件来选取方程中的参数，模拟得到的六边形斑图和条纹状斑图与实验照片很相似.

斑图动力学；Duffet-Boissonade 反应扩散方程组；数值模拟

斑图是在空间或时间上具有某种规律性的非均匀宏观结构.自然界的斑图可分为两类：一类是在热力学平衡态条件下的斑图，如无机化学中的晶体结构、有机聚合物中自组织形成的斑图，它的形成机理可以用平衡态热力学和统计物理原理来解释；另一类是在偏离热力学平衡态的条件下产生的斑图，如天上的条状云、水面上的波浪、动物体皮肤表面的花纹等，它们不能用平衡态热力学原理来解释，但可以从动力学角度来进行研究.斑图动力学就是研究这类斑图形成的原因及其规律的学科.

1952年英国数学家图灵在论文《形态形成的化学基础》[1]中，首次提出一个数学模型来描述生物体表面图案(如斑马身上的斑图)的生成过程.他设想，在生物体内存在一种“形态子”(比如是某种大分子)，在生物发育过程中，“形态子”与其他反应物发生了生物化学反应并在体内扩散.当 “形态子”浓度在空间分布变得不均匀时就会形成斑图.扩散本是一种常见的物理现象，如不同浓度物质之间的扩散，它们并不会形成斑图.为何图灵的反应扩散模型能产生斑图呢？因为他提出了一种新的假设：在扩散过程中存在两种相反的作用子，“活化子”与“阻滞子”，活化子催化反应，阻滞子阻碍反应，二者的空间传播速度不同，互相影响，造成空间各个区域的活化子与阻滞子浓度分布不同，影响到形态子的浓度在空间的分布也不同，于是形成了相对稳定的斑图.这种由体系内在的反应扩散特性所引起的空间均匀态被破坏的过程称为图灵失稳(图灵分岔)，所形成的图纹称为图灵斑图.

图灵斑图的研究直到上世纪60 年代末才开始引起重视，当时诺贝尔化学奖获得者普里高金(L.Prigogine) 为首的比利时布鲁塞尔热力学小组从热力学角度研究图灵斑图[6].他们证明，在远离热力学平衡态的条件下，体系可以通过自组织行为形成斑图，并称之为耗散结构，这是自然界不同系统中斑图形成的共性.同时，他们还提出了一个简单的、不违反任何化学反应动力学常识的反应模型——布鲁塞尔子(brusselator).此后科学家们对这个模型从理论、实验以及计算机模拟方面进行了大量的研究，大大推动了耗散理论以及斑图动力学的发展.

利用耗散结构理论研究化学振荡反应的另一个著名的例子是B-Z反应.B-Z反应是指1951年苏联生物学家Belousov 发现的一种在溶液中颜色会出现周期变化的化学反应，由于这个现象是违背热力学第二定律的，所以不被当时的科学界所接受.到了1961年研究生Zhabotinsky对该反应体系作了进一步的研究与改进，使溶液的颜色在红色和蓝色之间来回显著地变化，由于当时非线性理论已经有了突破，该反应体系也就逐步被人们所接受.该实验的发现是化学体系中的非单调行为概念的一个实验突破，它的诞生丰富了非线性动力学理论的内容，人们为了纪念该反应的发现者就把它称为B-Z反应.

在对布鲁塞尔子和B-Z反应等一系列模型的研究中，人们发现了大量的非线性现象.例如在时间序列上出现的多峰的周期振荡，迟豫振荡到阵发振荡的过渡，直到后来确认化学混沌现象和分岔现象的存在.在时空上也发现了多姿多彩的斑图，如靶形波，单螺旋波，双螺旋波，破碎螺旋波等[7].在《大学物理》杂志上也曾介绍过有关B-Z反应的实验[8].这些研究的进展使斑图动力学成了非线性动力学研究的重要方向.

为了更直观地了解斑图动力学，我们下面讨论一下《非线性科学与斑图动力学导论》[2]一书中介绍的一个简单的反应扩散模型即Duffet-Boissonade方程组，这是V.Duffet和J.J.Boissonade 在1992年发表在Chamical Physics上的论文《常规的与非常规的图灵斑图》[5]中首次建立的一个无量纲化的方程：

(1)

在方程(1)中u、v分别代表活化子与阻滞子的浓度，其中α、β、γ>0代表阻滞子与活化子的扩散系数比，d是阻滞子的扩散系数.

将这个方程组与扩散系数为a的物质扩散方程：

(2)

作个对比可知Duffet-Boissonade方程组中活化子的扩散系数是1，因此按照图灵斑图生成的条件，阻滞子的扩散系数d必须大于1.此外，它与物质扩散方程还有一个显著不同，就是在方程组中出现了u、v的交叉项，其中方程的三次项系数为负数，保证系统变量在失稳后不会因微扰而无限增长．二阶项的存在除去了方程在(u，v)→(-u，-v)变换下的不变性，这种对称性在化学系统中—般不存在；当γ=0时，可以回到这种对称的情形.这个模型的优点是可以对系统的线性行为与非线性行为分别处理．系统的均匀定态解的稳定条件由α、β、γ决定.在失稳后系统的行为由非线性项控制.

根据文献[2]中的介绍，方程有3个均匀定态解，当γ很小时，如果有1<β<α，系统在定态解u0=0，v0=0附近可能出现图灵分岔，要出现图灵分岔还要满足：

(d-β)2-4d(α-β)=0，d>β

当α<(β+d)2/4d时，图灵斑图开始生长.

文献[2]的分析对定性理解Duffet-Boissonade方程组很有帮助，但是由于没有给出这个方程组的数值模拟，很难想象方程组所呈现的复杂图像.本文利用文献[2]给出的图灵斑图的生成条件，用MATLAB计算了方程组在一维与二维情形的数值解，并将结果可视化，显示了这种反应扩散方程描述的斑图出现的过程与图案，这些斑图与实验出现的某些图像有一定的相似性.我们的模拟计算表明，文献[2]中对控制参数的分析基本正确，必须根据图灵分岔的条件去选择参数才能生成斑图，另外我们发现，初始条件决定了斑图的具体形状，不同的初始条件会生成不同的斑图，斑图稳定以后能维持很长时间，一般不会转化成别的斑图类型.这些结果证实Duffet-Boissonade模型确实能够描述反应扩散中的基本现象.

斑图动力学是非线性物理的重要分支，适当地了解斑图动力学可以帮助学生开阔视野，提高兴趣.本文的数值模拟所需要的知识仅仅是非线性偏微分方程组，对于学过计算物理学的学生而言，并不难掌握.这种学习既能帮助形象地理解与学习斑图动力学，也是计算物理在教学科研中的一个应用.

本文的计算都假定在区域的边界上存在无通量边界条件[4]，即函数在边界上法向导数为零，表示物质没有流出与流入.初始条件则是在区域存在一个小的微扰，微扰的传播形成图灵斑图.在将偏微分方程离散化时，对时间变量的一阶偏微分取向前差分公式，对空间变量的二阶差分取中心差分公式[3].计算中的参数按照前面介绍的图灵分岔的条件取值，为了对比，二维与一维的计算取相同的参数，一律取α=1.12，β=1.1，γ=0.01，d=1.5.

1 一维情形

计算区域为0≤x≤100，x的步长为1，时间t步长也取1，为了表现斑图生成的过程，步数t取了从小到大的6个值，分别为t=1， 15， 100，200，1000，15000.计算表明，步数超过15000以后，斑图形状基本稳定.初始条件为区域中心点存在一个小微扰.为了对比我们将u(实线)和v(点线)的图形画在同一张图中.

可以看到t=15时，原来为正值的微扰点两侧各出现了一个较小的负值的微扰；随着步数t增加，振荡以波动方式传开，振幅逐渐增大，振幅量级由最初的10-6增大到10-1；当步数t足够大时(例如到15 000以后)，振幅趋于平稳，形成稳定的条纹.图中u、v的波型基本重合，只是v波的振幅都要小于u波.这表明在相同位置上的活化子浓度小于阻滞子浓度，这正好符合图灵给出的生成图灵斑图的假设.

图1u，v的浓度分布曲线.

2 二维情形

文献[2]中登载有几张实验图片如图2所示，其中标志性的特征是存在六边形的点状斑图与直线或曲线型的条纹斑图.上面计算的一维情形虽然简单，但是不便与实验对比.为此，下面再作二维的数值模拟.

图2 实验中获取的斑图照片

下面计算中二维平面区域取为80×80.取时间步数t=10，200，1 000，15 000.计算考虑了4种不同的初始条件：1) 即中心点的小扰动，2) 在整个平面区域中均匀分布的点扰动，3) 直线形的扰动，4) 平行线状的扰动.由于u、v的图形很相似，所以只画出了u的图形.

1) 初始条件为平面区域的中心点有一小扰动.

开始出现近似于同心圆的斑图，随后图形轮廓越来越像六边形.这个结果对应于实验上所观测到的六边形图案.

2) 初始条件是在区域中均匀分布的点状的小扰动.

在斑图形成的过程中，初始扰动从各个扰动点中心对称地向外扩散，初期整个区域仍然有点状的分布，经过足够长的时间后有些点联成了条纹.注意它们始终不会成为均匀的分布.这些点状分布的图形4(a)、(b)与实验图片2(a)、(c)相似.

3) 初始的扰动是一条直线.

经过扩散形成的斑图是一系列平行线状条纹，这种图案形成后便保持基本稳定.它与实验图片2(b)有些相似.

4) 初始扰动是随机的浓度分布.

这个结果究竟如何很难想象.计算模拟表明，开始是呈现不规则的点状的分布，经过足够长的时间，也能形成条纹，但是其形状无法预测.下图显示的条纹是弯曲的，与实验图片2(d)有些相似.在数值模拟中也出现过比较直的条纹.

图3 图(a)、(b)、(c)、(d)是初始条件为区域中心点有一点状扰动在步数t=10，200，1000，15000时所生成的u的图形

图4 图(a)、(b)、(c)、(d)是初始条件为区域中心有均匀分布的点状扰动在时间步数t=10，200，1000，15000时所生成的u的图像

图5 图(a)、(b)、(c)、(d)是初始条件为区域中部有一直线形扰动时在时间步数t=10，200，1000，15000时所生成的u的图像

图6 图(a)、(b)、(c)、(d)是初始条件为区域内随机分布浓度在步数t=10，200，1000，15000时所生成的u的图像

3 讨论

通过数值模拟，我们发现有以下几点值得注意.

1) Duffet-Boissonade方程组虽然简单也能比较好地描述反应扩散现象的基本特性，按照稳定性分析的条件选取参数就能模拟出丰富的反应扩散现象，如同我们前面所展示的那样.但是，我们在数值计算中也发现文献[2]给出条件(d-β)2-4d(α-β)=0 并不严格.根据这个条件和我们选取的α、β值，似乎d的取值只能是1.439 33，而我们的数值模拟表明，d的取值从1.3到1.8都是可行的.

2) 文献[2]在解释活化子与阻滞子对斑图形成的作用机制也与我们数值模拟的结果不一致.该文献使用了如下图形描述活化子与阻滞子的传播过程.

图7 文献[2]中表示活化子与阻滞子扩散速度的示意图

文献[2]认为由于阻滞子的扩散速度远大于活化子的，阻滞子将会更快地经扩散离开微扰形成上述图形，注意这里有间断的波形.但在我们的数值模拟中，无论u还是v的波形，在传播中都是连续的，不会发生间断.更重要的是二者的波形基本相似，只是振幅不同.为此我们设想了斑图形成机制的另一种解释：在传播过程中，二者的波形其实是相同的(见图1)，但由于振幅不同造成二者在各点的浓度不同.活化子是催化反应，而阻滞子是减缓反应.二者的浓度差决定了在各处的反应生成物的浓度也不同，浓度不同就会产生斑图.

3) 初始条件会决定斑图的形状.不同初始条件会生成不同斑图，不改变控制参数则它们不会互相转化.这个结果给我们印象很深.

斑图动力学仍然是目前非线性物理中热门的研究方向，有许多模型在进行理论、实验和数值模拟的研究.本文的数值模拟提供了一些有趣的初步结果，限于所选取的符合分岔条件的参数值非常有限，所以还有许多工作值得深入去做，这可以作为今后研究的内容.

[1] Turing A M.The Chemical basis of morphogenesis[J]. Phil， Trans. R. Soc. London Ser B， 327， 37(1952).

[2] 欧阳颀.非线性科学与斑图动力学导论[M].北京：北京大学出版社，2010:156，138-164.

[3] 彭芳麟.计算物理基础[M].北京：高等教育出版社，2010:330-344.

[4] 刘迎东.图灵斑图动力学的数学机制[J].北方交通大学学报，2004，28(3):1-3.

[5] Duffet V， Boissonade J J. Conventional and unconventional Turing patterns[J]. Chem Phys，1992，96:664.

[6] Nicolis G ，Prigongine 1. Self- Organizalion in Nonequilibriwn Chemical Syslems[M] . New York: Wiley,1977.

[7] 宗春燕，王玉梅，高庆宇.Belousov-Zhabotinsky反应研究进展[J].淮阴工学院学报， 2006，15(1):54-57.

[8] 郑佳喻，周华喜，孙金芳，等. B-Z反应中斑图的数值模拟及实验研究[J].大学物理，2008, 27 (10)：50-54.

Numerical simulation of the Duffet-Boissonade equations in pattern dynamics

PENG Fang-lin，JU Xiao-lin, LIU Chang

(Department of Physics, Beijing Normal University, Beijing 100875, China)

The one- and two-dimensional numerical simulations of Duffet-Boissonade reaction-diffusion equations in pattern dynamics are performed. Setting up four kinds of different initial conditions and selecting the parameters in the equation according to Turing bifurcation condition, the hexagon pattern and fringe pattern obtained from the simulation are very similar to the experimental results.

pattern dynamics；Duffet-Boissonade reaction-diffusion equations ; numerical simulation

2016-03-28；

2016-05-26

彭芳麟(1947—)，男，江西泰和人，北京师范大学物理学系教授，主要从事计算物理教学和研究工作.

教学研究

O 411.3

A

1000- 0712(2016)12- 0001- 06