参数依赖延迟的Gray⁃Scott模型稳定性和分支分析

Stability and Hopf Bifurcations for Gray-Scott Model with Delay-Dependent Parameters

摘要

研究具有延迟反馈控制的Gray⁃Scott模型稳定性和分支问题，其中反馈控制的强度是延迟\tau的非线性函数{\mathrm{e}}^{-p\tau }。选择(p,\tau )为参数，结合几何方法和解析方法，确定系统稳定性转换和Hopf分支的临界值。随着两个参数变化，得到Hopf分支曲线和稳定区域。最后数值模拟验证所得结论。

Abstract

The stability and bifurcation of a Gray-Scott model with delay feedback control are investigated， in which the strength of the feedback control is decided by a nonlinear function of the delay τ， e^-pτ . By choosing decay rate p and time delay τ as the parameters combined with the geometric and analytical methods， the critical values of the system stability conversion and Hopf bifurcation are determined respectively. The Hopf bifurcation curve and the stability region are obtained by changing two parameters. Finally， the conclusions are confirmed through numerical simulation.

Graphical abstract

关键词

Gray⁃Scott模型 / 延迟 / 稳定性转换 / Hopf分支

Key words

Gray-Scott model / delay / stability conversion / Hopf bifurcation

Gray⁃Scott模型^［1⁃2］是一种重要的反应扩散系统，起源于模拟自催化化学反应，其化学反应过程为

\begin{array}{l}\mathrm{U}+\mathrm{2}V\to \mathrm{3}V,\\ V\to P,\end{array}

其中U,V为两种化学物质，V作为催化剂对自身生产起催化作用，P为一种惰性产物，不参与反应过程。该反应过程的无量纲反应扩散方程形式为

\left\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}=D_u{\nabla }^{\mathrm{2}}u+a(\mathrm{1}-u)-u{v}^{\mathrm{2}},\\ \frac{\partial v(x,t)}{\partial t}=D_v{\nabla }^{\mathrm{2}}v-(a+b)v+u{v}^{\mathrm{2}},\end{array}\right.（1）

其中u(x,t),v(x,t)分别表示在x点处t\ge \mathrm{0}时抑制剂与催化剂的浓度，a表示入流率，a+b表示催化剂V的移出率，D_u和D_v分别表示化学物质U,V的扩散系数。

迄今，已有许多文章从不同的角度对系统（1）进行了研究^［3⁃12］。2009年，Kyrychko等^［11］对系统（1）施加了延迟反馈控制项

\frac{\partial }{\partial t}\left(\begin{array}{l}u\\ v\end{array}\right)=\left(\begin{array}{l}f(u,v)\\ g(u,v)\end{array}\right)+\left(\begin{array}{cc}{D}_u& \mathrm{0}\\ \mathrm{0}& D_v\end{array}\right){\nabla }^{\mathrm{2}}\left(\begin{array}{l}u\\ v\end{array}\right)+KA\left(\begin{array}{l}u(t-\tau )-u(t)\\ v(t-\tau )-v(t)\end{array}\right),（2）

其中f(u,v)=-u{v}^{\mathrm{2}}+a(\mathrm{1}-u),g(u,v)=u{v}^{\mathrm{2}}-(a+b)v，K为反馈强度，A为2×2的控制矩阵。研究控制矩阵A取\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right)、\left(\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right)和\left(\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\varphi & \mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ -\mathrm{s}\mathrm{i}\mathrm{n}\varphi & \mathrm{c}\mathrm{o}\mathrm{s}\varphi \end{array}\right)(\varphi \mathrm{是相}\mathrm{位})等情况时，延迟反馈控制对系统（2）的时空斑图的影响。在每种情况下，在控制强度和时滞的参数空间中找到了稳定边界。在时空混沌的情况下，施加的控制既可以稳定均匀稳态，也可以导致平凡稳态和传播行波之间的双稳态。

此外，为了获得关于系统（1）的动力学的信息，研究者也对常微分系统进行了讨论。

\left\{\begin{array}{c}\frac{\mathrm{d}u}{\mathrm{d}t}=-u{v}^{\mathrm{2}}+a(\mathrm{1}-u),\\ \frac{\mathrm{d}v}{\mathrm{d}t}=-(a+b)v+u{v}^{\mathrm{2}}\mathrm{。}\end{array}\right.（3）

系统（3）是系统（1）中D_u=D_v=\mathrm{0}的情形。文献［12］研究了系统（3）对应的三次多项式系统在庞加莱盘中的相位图。相位图和相应的分岔图显示了这类系统的丰富性和复杂的动力学。

基于上述讨论，本文将研究如下形式的一类具有延迟反馈项Gray⁃Scott模型：

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}t}=-u{v}^{\mathrm{2}}+a(\mathrm{1}-u)+K{\mathrm{e}}^{-p\tau }(v(t-\tau )-v(t)),t>\mathrm{0},\\ \frac{\mathrm{d}v}{\mathrm{d}t}=-(a+b)v+u{v}^{\mathrm{2}},t>\mathrm{0},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ u(t)=\varphi (t),\hspace{1em}v(t)=\phi (t),t\in [-\tau ,\mathrm{0}],\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.（4）

其中，p>\mathrm{0},\tau >\mathrm{0}表示延迟；K是常数，K{\mathrm{e}}^{-p\tau }表示延迟反馈强度。受过去历史记忆消退效应的启发，延迟反馈强度随着延迟\tau的增加而衰减，当延迟\tau趋于无穷时，反馈就消失。对于依赖于延迟的反馈控制研究也已取得了很多成果^{［13⁃16］}。通过研究成果可以看出，参数对延迟\tau的依赖性使得系统的动力学性质变的更加复杂。本文选择延迟\tau和衰减率p作为分岔参数，将分析系统（4）产生Hopf分岔和稳定转换的条件。

1 正平衡点的稳定性和Hopf分支的存在性

经过计算，系统（4）存在三个平衡点，分别为

\begin{array}{l}{E}_{\mathrm{0}}=(u_{\mathrm{0}},v_{\mathrm{0}})=(\mathrm{1,0}),\\ E_{\mathrm{1}}=(u_{\mathrm{1}},v_{\mathrm{1}})=(\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}-\sqrt[]{A}),\frac{\mathrm{1}}{\mathrm{2}}\cdot \frac{a}{a+b}\cdot (\mathrm{1}+\sqrt[]{A})),\\ E_{\mathrm{2}}=(u_{\mathrm{2}},v_{\mathrm{2}})=(\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}+\sqrt[]{A}),\frac{\mathrm{1}}{\mathrm{2}}\cdot \frac{a}{a+b}\cdot (\mathrm{1}-\sqrt[]{A})),\end{array}

其中A=\mathrm{1}-\frac{\mathrm{4}{(a+b)}^{\mathrm{2}}}{a}\mathrm{。}根据文献［11］可知常值平衡点E_{\mathrm{0}}总是稳定的，而平衡点E_{\mathrm{2}}对于它存在的任何参数的齐次扰动都是不稳定的，因此本文选取平衡点E_{\mathrm{1}}进行研究。

在平衡点E_{\mathrm{1}}处将系统（4）线性化，得

\left(\begin{array}{l}\frac{\partial u(x,t)}{\partial t}\\ \frac{\partial v(x,t)}{\partial t}\end{array}\right)=L_{\mathrm{1}}\left(\begin{array}{l}u(x,t)\\ v(x,t)\end{array}\right)+L_{\mathrm{2}}\left(\begin{array}{l}u(x,t-\tau )\\ v(x,t-\tau )\end{array}\right),

L_{\mathrm{1}}=\left(\begin{array}{ll}-a-v_{\mathrm{1}}^{\mathrm{2}}& -(\mathrm{2}(a+b)+K{\mathrm{e}}^{-p\tau })\\ v_{\mathrm{1}}^{\mathrm{2}}& a+b\end{array}\right),L_{\mathrm{2}}=\left(\begin{array}{ll}\mathrm{0}& K{\mathrm{e}}^{-p\tau }\\ \mathrm{0}& \mathrm{0}\end{array}\right)\mathrm{。}

系统（4）的特征方程为

{\lambda }^{\mathrm{2}}+(v_{\mathrm{1}}^{\mathrm{2}}-b)\lambda +(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }-K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\cdot {\mathrm{e}}^{-\lambda \tau }=\mathrm{0}\mathrm{。}（5）

令\lambda =\mathrm{i}\omega，分离实虚部，得到

\left\{\begin{array}{l}-{\omega }^{\mathrm{2}}+(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }-K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )=\mathrm{0},\\ (v_{\mathrm{1}}^{\mathrm{2}}-b)\omega +K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )=\mathrm{0}\mathrm{。}\end{array}\right.（6）

解得

\left\{\begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )=\frac{-{\omega }^{\mathrm{2}}+(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }}{K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }}=\frac{-{\omega }^{\mathrm{2}}+(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)}{K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }}+\mathrm{1},\\ \mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )=\frac{-(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega }{K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }}\mathrm{。}\end{array}\right.（7）

令\theta =\omega \tau，由于\tau \in [\mathrm{0},+\mathrm{\infty })，从而\theta \in [\mathrm{0},+\mathrm{\infty })\mathrm{。}由式（7）可推得

\frac{\mathrm{c}\mathrm{o}\mathrm{s}\theta -\mathrm{1}}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }=\frac{-{\omega }^{\mathrm{2}}+(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)}{-(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega }\mathrm{。}（8）

化简得

\mathrm{s}\mathrm{i}\mathrm{n}\theta \cdot {\omega }^{\mathrm{2}}+(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\theta )(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega -\mathrm{s}\mathrm{i}\mathrm{n}\theta (a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)=\mathrm{0}\mathrm{。}（9）

解得

{\omega }_{\pm }=\frac{-(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\theta )(v_{\mathrm{1}}^{\mathrm{2}}-b)\pm \sqrt[]{[(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\theta )(v_{\mathrm{1}}^{\mathrm{2}}{-b)]}^{\mathrm{2}}+\mathrm{4}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta (v_{\mathrm{1}}^{\mathrm{2}}-a)(a+b)}}{\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\theta }\mathrm{。}（10）

显然，若，有\mathrm{s}\mathrm{i}\mathrm{n}\theta >\mathrm{0},{\omega }_{+}>\mathrm{0}；若K>\mathrm{0},v_{\mathrm{1}}^{\mathrm{2}}>\mathrm{m}\mathrm{a}\mathrm{x}(a,b)，有。本文仅考虑情况（K>\mathrm{0},v_{\mathrm{1}}^{\mathrm{2}}>\mathrm{m}\mathrm{a}\mathrm{x}(a,b)讨论方法类似）。

根据式（7），有

{\mathrm{e}}^{-p_{+}\tau }=\frac{-(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega }{K{v}_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )}\mathrm{。}（11）

进一步，有

\left\{\begin{array}{l}{p}_{+}=-\frac{{\omega }_{+}}{\theta }\mathrm{l}\mathrm{n}(\frac{-(v_{\mathrm{1}}^{\mathrm{2}}-b){\omega }_{+}}{K{v}_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\theta }),\\ {\tau }_{+}=-\frac{\mathrm{1}}{p_{+}}\mathrm{l}\mathrm{n}(\frac{-(v_{\mathrm{1}}^{\mathrm{2}}-b){\omega }_{+}}{K{v}_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\theta })\mathrm{。}\end{array}\right.（12）

令\theta ={\theta }_{\mathrm{0}}+\mathrm{2}n\mathrm{\pi },n=\mathrm{0,1},\mathrm{2}\cdots，由于\mathrm{s}\mathrm{i}\mathrm{n}\theta >\mathrm{0}，那么{\theta }_{\mathrm{0}}\in (\mathrm{0},\mathrm{\pi })。

当\tau =\mathrm{0}时，方程（5）变为

{\lambda }^{\mathrm{2}}+(v_{\mathrm{1}}^{\mathrm{2}}-b)\lambda +(a+b)(v_{\mathrm{1}}^{\mathrm{2}}-a)=\mathrm{0}\mathrm{。}

根据韦达定理有，可知系统在\tau =\mathrm{0}时是稳定的。

讨论\tau >\mathrm{0}的情况。由于\theta满足式（12），且减小p的值或增大\tau的值会导致\theta的个数持续增加，这种变化同时也会引起Hopf分支。由式（12）中p=p_{+}的几何图，可以直观看到Hopf分支的出现频率，如图1所示，其中n=\mathrm{0,1},\mathrm{2,3}的情况分别由蓝色、橙色、黄色、紫色线条表示。根据图1可得以下结论：

（1） 若p>p_{\mathrm{0}}^{\mathrm{m}\mathrm{a}\mathrm{x}}，特征方程（5）不存在纯虚根，此时不产生Hopf分支，系统（4）对所有的\tau >\mathrm{0}都是稳定的；

（2） 若，特征方程（5）存在\mathrm{2}(n+\mathrm{1})个纯虚根，这表明随着延迟的改变系统（4）会发生\mathrm{2}(n+\mathrm{1})次Hopf分支；

（3） 若p=p_n^{\mathrm{m}\mathrm{a}\mathrm{x}}，特征方程（5）存在（\mathrm{2}n+\mathrm{1}）个纯虚根，这表明随着时延迟的改变系统（4）会发生（\mathrm{2}n+\mathrm{1}）次Hopf分支。

对每个n，都能通过\frac{\mathrm{d}{p}_{+}(\theta )}{\mathrm{d}\theta }=\mathrm{0}计算最大值p_n^{\mathrm{m}\mathrm{a}\mathrm{x}}。由式（12），有

\frac{\mathrm{d}{p}_{+}}{\mathrm{d}\theta }=-\frac{p_{+}}{\theta }+\frac{p_{+}\tau }{\theta }\cdot \frac{\mathrm{d}\omega }{\mathrm{d}\theta }-\frac{\mathrm{1}}{\theta }\cdot \frac{\mathrm{d}\omega }{\mathrm{d}\theta }+\frac{\omega }{\theta \mathrm{t}\mathrm{a}\mathrm{n}\theta }\mathrm{。}（13）

而由式（9），有

\frac{\mathrm{d}\omega }{\mathrm{d}\theta }=\frac{-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{s}\mathrm{e}{\mathrm{c}}^{\mathrm{2}}\frac{\theta }{\mathrm{2}}\cdot (v_{\mathrm{1}}^{\mathrm{2}}-b)\cdot \omega }{\mathrm{2}\omega +\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}\cdot (v_{\mathrm{1}}^{\mathrm{2}}-b)}\mathrm{。}（14）

结合式（13）与式（14），有

\frac{\mathrm{d}{p}_{+}}{\mathrm{d}\theta }=\frac{M}{\mathrm{2}\theta \mathrm{t}\mathrm{a}\mathrm{n}\theta (\mathrm{2}\omega +\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b))},（15）

其中

\begin{array}{l}M=-\mathrm{2}p\mathrm{t}\mathrm{a}\mathrm{n}\theta (\mathrm{2}\omega +(v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}})-(p\tau -\mathrm{1})\mathrm{t}\mathrm{a}\mathrm{n}\theta \mathrm{s}\mathrm{e}{\mathrm{c}}^{\mathrm{2}}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega +\\ \mathrm{2}\omega (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}+\mathrm{4}{\omega }^{\mathrm{2}}\mathrm{。}\end{array}

将式（9）与式（10）代入式（15）中，并且运用\frac{\mathrm{d}p(\theta )}{\mathrm{d}\theta }=\mathrm{0}解出{\theta }_n^{\mathrm{m}\mathrm{a}\mathrm{x}}，通过式（12）可得到p_n^{\mathrm{m}\mathrm{a}\mathrm{x}}的值。图1给出n取不同值时({\theta }_n^{\mathrm{m}\mathrm{a}\mathrm{x}},p_n^{\mathrm{m}\mathrm{a}\mathrm{x}})的位置。此外，时滞\tau的Hopf分支值也可以通过式（12）求得。对于给定的p，将方程p(\theta )=p的解定义为{\theta }_n，其对应的特征根记作{\omega }_n，对应的时滞值记作{\tau }_n。

讨论系统（4）的稳定性转换和Hopf分支。对于给定的\tau，假设方程（5）存在一组共轭根\lambda =\alpha (\tau )\pm \mathrm{i}\omega (\tau )。为了简单起见，定义{\tau }_c={\tau }_n,{\omega }_c={\omega }_n(n=\mathrm{0,1},\mathrm{2}\cdots )，此时有\alpha ({\tau }_c)=\mathrm{0}以及\omega ({\tau }_c)={\omega }_c，对方程（5）关于\tau求导，有

\begin{array}{l}((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}\tau {\mathrm{e}}^{-p\tau }\mathrm{c}\mathrm{o}\mathrm{s}\theta )\frac{\mathrm{d}\alpha }{\mathrm{d}\tau }+(K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\tau \mathrm{s}\mathrm{i}\mathrm{n}\theta -\mathrm{2}{\omega }_n)\frac{\mathrm{d}{\omega }_n}{\mathrm{d}\tau }=pK{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }-K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }(p\mathrm{c}\mathrm{o}\mathrm{s}\theta +{\omega }_n\mathrm{s}\mathrm{i}\mathrm{n}\theta ),\\ (\mathrm{2}{\omega }_n-\tau K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\mathrm{s}\mathrm{i}\mathrm{n}\theta )\frac{\mathrm{d}\alpha }{\mathrm{d}\tau }+((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}\tau {\mathrm{e}}^{-p\tau }\mathrm{c}\mathrm{o}\mathrm{s}\theta )\frac{\mathrm{d}{\omega }_n}{\mathrm{d}\tau }=pK{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\mathrm{s}\mathrm{i}\mathrm{n}\theta -{\omega }_{n}K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p\tau }\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{。}\end{array}解得

\frac{\mathrm{d}\alpha }{\mathrm{d}\tau }{|}_{\tau ={\tau }_c}=\frac{F}{((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c{\mathrm{e}}^{-p{\tau }_c}{\mathrm{c}\mathrm{o}\mathrm{s}\theta )}^{\mathrm{2}}+(\mathrm{2}{\omega }_c-{\tau }_{c}K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}{\mathrm{s}\mathrm{i}\mathrm{n}\theta )}^{\mathrm{2}}},（16）

其中

\begin{array}{l}F=((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c{\mathrm{e}}^{-p{\tau }_c}\mathrm{c}\mathrm{o}\mathrm{s}\theta )(pK{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}-K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}(p\mathrm{c}\mathrm{o}\mathrm{s}\theta +{\omega }_c\mathrm{s}\mathrm{i}\mathrm{n}\theta )-\\ (K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}{\tau }_c\mathrm{s}\mathrm{i}\mathrm{n}\theta -\mathrm{2}{\omega }_c)(pK{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}\mathrm{s}\mathrm{i}\mathrm{n}\theta -{\omega }_{c}K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}\mathrm{c}\mathrm{o}\mathrm{s}\theta )=\\ {\mathrm{e}}^{-\mathrm{2}p{\tau }_c}\cdot K{v}_{\mathrm{1}}^{\mathrm{2}}((v_{\mathrm{1}}^{\mathrm{2}}-b){\mathrm{e}}^{p{\tau }_c}+K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c\mathrm{c}\mathrm{o}\mathrm{s}\theta )(p-p\mathrm{c}\mathrm{o}\mathrm{s}\theta -{\omega }_c\mathrm{s}\mathrm{i}\mathrm{n}\theta )-\\ K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-\mathrm{2}p{\tau }_c}(K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c\mathrm{s}\mathrm{i}\mathrm{n}\theta -\mathrm{2}{\omega }_c{\mathrm{e}}^{p{\tau }_c})(p\mathrm{s}\mathrm{i}\mathrm{n}\theta -{\omega }_c\mathrm{s}\mathrm{i}\mathrm{n}\theta )\mathrm{。}\end{array}（17）

将式（11）代入式（16），化简得

\begin{array}{l}F=\frac{{\mathrm{e}}^{-\mathrm{2}p{\tau }_c}(K{v}_{\mathrm{1}}^{\mathrm{2}}{)}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta }{\mathrm{2}{\omega }_c(v_{\mathrm{1}}^{\mathrm{2}}-b)}(-\mathrm{2}p\mathrm{t}\mathrm{a}\mathrm{n}\theta \cdot ((v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}+\mathrm{2}{\omega }_c)-\\ \mathrm{2}p{\omega }_c{\tau }_c(v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{c}\mathrm{o}\mathrm{s}\theta }+\mathrm{2}{\omega }_c(v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\theta +\mathrm{4}{\omega }_c^{\mathrm{2}})\mathrm{。}\end{array}

同时，由于

\begin{array}{l}-(p\tau -\mathrm{1})\mathrm{t}\mathrm{a}\mathrm{n}\theta \mathrm{s}\mathrm{e}{\mathrm{c}}^{\mathrm{2}}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega +\mathrm{2}\omega (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}=\\ \mathrm{2}p\tau \omega (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{c}\mathrm{o}\mathrm{s}\theta }+\mathrm{2}(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega \mathrm{t}\mathrm{a}\mathrm{n}\theta ,\end{array}

故

\begin{array}{l}-\mathrm{2}p\mathrm{t}\mathrm{a}\mathrm{n}\theta \cdot ((v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}+\mathrm{2}\omega )-\mathrm{2}p\omega \tau (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{c}\mathrm{o}\mathrm{s}\theta }+\mathrm{2}\omega (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\theta +\mathrm{4}{\omega }^{\mathrm{2}}=\\ -\mathrm{2}p\mathrm{t}\mathrm{a}\mathrm{n}\theta (\mathrm{2}\omega +(v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}})-(p\tau -\mathrm{1})\mathrm{t}\mathrm{a}\mathrm{n}\theta \mathrm{s}\mathrm{e}{\mathrm{c}}^{\mathrm{2}}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b)\omega +\\ \mathrm{2}\omega (v_{\mathrm{1}}^{\mathrm{2}}-b)\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}+\mathrm{4}{\omega }^{\mathrm{2}}=M\mathrm{。}\end{array}

则

\begin{array}{l}\frac{\mathrm{d}\alpha }{\mathrm{d}\tau }\cdot \frac{\mathrm{d}p}{\mathrm{d}\theta }{|}_{\tau ={\tau }_c}=\frac{{\mathrm{e}}^{-\mathrm{2}p{\tau }_c}(K{v}_{\mathrm{1}}^{\mathrm{2}}{)}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta }{\mathrm{2}{\omega }_c(v_{\mathrm{1}}^{\mathrm{2}}-b)}\cdot \frac{M}{((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c{\mathrm{e}}^{-p{\tau }_c}{\mathrm{c}\mathrm{o}\mathrm{s}\theta )}^{\mathrm{2}}+(\mathrm{2}{\omega }_c-{\tau }_{c}K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}{\mathrm{s}\mathrm{i}\mathrm{n}\theta )}^{\mathrm{2}}}\cdot \\ \frac{M}{\mathrm{2}\theta \mathrm{t}\mathrm{a}\mathrm{n}\theta (\mathrm{2}{\omega }_c+\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b))}=\frac{(K{v}_{\mathrm{1}}^{\mathrm{2}}{)}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta }{\mathrm{4}\theta {\omega }_c(v_{\mathrm{1}}^{\mathrm{2}}-b)(\mathrm{2}{\omega }_c+\mathrm{t}\mathrm{a}\mathrm{n}\frac{\theta }{\mathrm{2}}(v_{\mathrm{1}}^{\mathrm{2}}-b))}\cdot \\ \frac{M^{\mathrm{2}}}{((v_{\mathrm{1}}^{\mathrm{2}}-b)+K{v}_{\mathrm{1}}^{\mathrm{2}}{\tau }_c{\mathrm{e}}^{-p{\tau }_c}{\mathrm{c}\mathrm{o}\mathrm{s}\theta )}^{\mathrm{2}}+(\mathrm{2}{\omega }_c-{\tau }_{c}K{v}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{e}}^{-p{\tau }_c}{\mathrm{s}\mathrm{i}\mathrm{n}\theta )}^{\mathrm{2}}}>\mathrm{0}\mathrm{。}\end{array}

由图1可知，对于每一个n，p_n^{\mathrm{m}\mathrm{a}\mathrm{x}}都是唯一的，因此，若\frac{\mathrm{d}p}{\mathrm{d}\theta }{|}_{\theta ={\theta }_c}>\mathrm{0},则有\frac{\mathrm{d}\alpha }{\mathrm{d}\tau }{|}_{\tau ={\tau }_c}>\mathrm{0};若则有

结合之前的研究，可知系统（4）在\tau =\mathrm{0}时是稳定的。当\tau继续增大时，情况保持不变，直到对应一对实部为零的特征根的第一个临界值{\tau }_{\mathrm{0}}的出现，此时系统开始由稳定变为不稳定状态，并在{\tau }_{\mathrm{0}}处产生Hopf分支；当\tau继续增大到第二个临界值{\tau }_{\mathrm{1}}时，系统由不稳定再次变为稳定状态。因此在(\tau ,p)平面内形成了一个稳定域，如图2所示，阴影部分即为稳定域，并且在这个区域的边界上会产生Hopf分支。

2 数值模拟

令a=\mathrm{0.1},b=\mathrm{0.02},K=-\mathrm{0.5},p=\mathrm{0.1}，此时平衡点E_{\mathrm{1}}=(u_{\mathrm{1}},v_{\mathrm{1}})=(\mathrm{0.174,0.688})。计算可知满足假设条件\mathrm{s}\mathrm{i}\mathrm{n}\theta >\mathrm{0},v_{\mathrm{1}}^{\mathrm{2}}>\mathrm{m}\mathrm{a}\mathrm{x}(a,b)。基于前文研究，当\tau =\mathrm{1}时，E_{\mathrm{1}}是稳定的，如图3所示；当\tau =\mathrm{2.514}时，产生周期解，如图4所示。此外，随着\tau的继续增大，解将趋于常值平衡点E_{\mathrm{0}}。另外，根据文献［11］可知在远离图灵和Hopf区域的参数范围内，系统总是以E_{\mathrm{0}}或E_{\mathrm{1}}结束。在图5给出了\tau =\mathrm{3}时，系统在E_{\mathrm{0}}处结束。

3 小结

本文研究了一类具有延迟反馈控制的Gray⁃Scott模型，其中反馈控制强度为K{\mathrm{e}}^{-p\tau }。这种参数依赖于延迟的反馈控制，使得系统的稳定性以及分支分析变得更加复杂。研究给出了确定平衡点稳定条件的几何方法和解析方法，并讨论了该模型的Hopf分支。主要结果如下：

（1） 除时滞\tau外，系统（4）的Hopf分支强烈依赖于参数p，虚特征根的横截条件由参数p(\theta )变化决定；

（2） 在(\tau ,p)平面上作出了系统（4）的Hopf分支曲线，并且标记出了稳定域（图2）；

（3） 本文通过数值模拟的方法验证了当\tau持续增大时，系统（4）会由稳定状态转为不稳定状态。

参考文献

[1] GRAY P， SCOTT S K.Autocatalytic reactions in the isothermal， continuous stirred tank reactor： Isolas and other forms of multistability［J］.Chemical Engineering Science，1983，38（1）：29-43.

[2] SCOTT S K， BODDINGTON T， GRAY P.Analytical expressions for effectiveness factors in non-isothermal spherical catalysts［J］.Chemical Engineering Science，1984，39（6）：1087-1097.

[3] NISHIURA Y， UEYAMA D.Spatio-temporal chaos for the Gray⁃Scott model［J］.Physica D：Nonlinear Phenomena，2001，150（3）：137-162.

[4] HAN C， WANG Y， LI Z.A high-precision numerical approach to solving space fractional Gray-Scott model［J］.Applied Mathematics Letters，2022，125：107759.

[5] WANG J， JIA Y.Analysis on bifurcation and stability of a generalized Gray-Scott chemical reaction model［J］.Physica A： Statistical Mechanics and Its Applications，2019，528：121394.

[6] FRANCESCA B， ROBERTO B， ELEONORA B，et al.Mathematical pattern for parametric design：The case study of Gray-Scott cross diffusion model［J］.Procedia Computer Science，2023，217：756-764.

[7] CHAI L， LIU Y， LI H，et al.Fast TT-M fourth-order compact difference schemes for a two-dimensional space fractional Gray-Scott model［J］.Computers & Mathematics with Applications，2023，141：191-206.

[8] LIANG J， JIANG N， LIU C，et al.On a reversible Gray-Scott type system from energetic variational approach and its irreversible limit［J］.Journal of Differential Equations，2022，309：427-454.

[9] ZHAI S， WENG Z， ZHUANG Q，et al.An effective operator splitting method based on spectral deferred correction for the fractional Gray-Scott model［J］.Journal of Computational and Applied Mathematics，2023，425：114959.

[10] SINGH S， MITTAL R C， RAHMAN THOTTOLI S，et al.High-fidelity simulations for turing pattern formation in multi-dimensional Gray-Scott reaction-diffusion system［J］.Applied Mathematics and Computation，2023，452：128079.

[11] KYRYCHKO Y N， BLYUSS K B， HOGAN S J，et al.Control of spatiotemporal patterns in the Gray-Scott model［J］. Chaos，2009，19（4）：043126.

[12] CHEN T， LI S， LLIBRE J.Phase portraits and bifurcation diagram of the Gray-Scott model［J］.Journal of Mathematical Analysis and Applications，2021，496（2）：124840.

[13] MA S， LU Q， FENG Z.Double Hopf bifurcation for van der Pol-Duffing oscillator with parametric delay feedback control［J］.Journal of Mathematical Analysis and Applications，2008，38（2）：993-1007.

[14] MA S， LU Q.Double Hopf bifurcation for Stuart-Landau system with nonlinear delay feedback and delay-dependent parameters［J］.Advances in Complex Systems，2007，10（4）：423-448.

[15] 罗佳伟，徐旭.参数与时滞相关延迟振子的稳定性及Hopf分岔分析［J］.吉林大学学报（理学版），2018，56（6）：1337-1344.

[16] 范丽，史忠科，陈斯养.参数依赖时滞的Nicholson生态模型的稳定性和分支［J］.应用数学学报，2010，33（5）：824-839.

基金资助

中央高校基本科研业务费专项资金资助项目“具有年龄结构的传染病模型的定性分析”(2572018BC19)