空间分数阶Gray-Scott方程的数值算法

Numerical Algorithm for the Gray-Scott Equation of Spatial Fractional Order

摘要: 本文基于算子分裂方法，提出了求解分数阶Gray-Scott模型的一种高效数值逼近格式。首先采用算子分裂法将原问题分解为线性子问题和非线性子问题：线性子问题采用Crank-Nicolson(CN)格式结合二阶中心差分，建立整体二阶的数值计算格式；非线性子问题采用CN格式结合Rubin-Graves线性化技术，建立线性化求解格式；并给出算法的稳定性和收敛性分析。最后，通过数值算例验证了算法的有效性。

Abstract: In this paper, an efficient numerical approximation algorithm for solving the fractional-order Gray-Scott model is proposed based on the operator splitting method. Firstly, the operator splitting method is used to decompose the original problem into linear and nonlinear subproblems: the lin-ear subproblem adopts the Crank-Nicolson (CN) format combined with the second-order central difference to establish the overall second-order numerical computation format; the nonlinear sub-problem adopts the CN format combined with the Rubin-Graves linearization technique to establish the linearized solution format; and the stability and convergence analysis of the algorithm are given. Finally, the validity of the algorithm is verified by numerical examples.

1. 引言

Gray-Scott模型(简称GS方程)是反应–扩散模型的重要组成部分，在化学、生物学、物理学等领域有广泛的应用。反应–扩散所形成的各种模式在自然界中都可以找到，例如点–点的分裂，条带以及运移波等。GS模型是科研学者描述反应器中浓度的时空变化的一种化学动力模型。整数阶GS方程如下：

$\{\begin{array}{l}\frac{\partial u}{\partial t}={\mu }_u{\Delta }_u-u{v}^2+F\left(1-u\right)\\ \frac{\partial v}{\partial t}={\mu }_v{\Delta }_v+u{v}^2-\left(F+\kappa \right)v\end{array}$

传统GS模型是一个二阶整数阶偏微分方程，但现实生活中许多自然现象无法用整数阶微积分方程模拟。例如，在研究分子扩散、粘性流体运动及物质运输过程等问题时，人们发现传统的基于Fick扩散定律建立的整数阶对流扩散方程，并不能很好的模拟溶质迁移过程中所检测到的穿透曲线提前和拖尾现象，以及弥散系数随研究尺度增大而增大的现象，而分数阶对流扩散方程 [1] [2] [3] [4] 可以很好的模拟这种现象。本文研究的分数阶GS方程如下：

$\{\begin{array}{l}\frac{\partial u}{\partial t}=-{\mu }_u{\left(-\Delta \right)}^{\frac{a}{2}}u-u{\nu }^2+F\left(1-u\right),\left(x,t\right)\in \Omega \times J\\ \frac{\partial v}{\partial t}=-{\mu }_v{\left(-\Delta \right)}^{\frac{\beta }{2}}v+u{\nu }^2-\left(F+\kappa \right)v,\left(x,t\right)\in \Omega \times J\end{array}$ (1.1)

分数阶偏微分方程被广泛应用于工程和科学领域的复杂系统中。其分数阶导数能够很好的描述具有遗传性或记忆性、长距离依赖性的复杂环境，因而分数阶偏微分方程成为了模拟不规则扩散、污染物运移、随机动态系统、经济以及风险测评等复杂现象的强有力的数学工具。

近年来，在工程与科学领域，许多研究人员对分数阶GS模型进行了研究和分析，但是大量研究表明求分数阶GS模型的精确解是极为困难的，因此一些学者转而研究数值计算技术求解GS模型。2019年，Wang等在文 [5] 中采用谱配置方法求解分数阶GS模型，并给出了相应的数值模拟。在文 [6] 中，Abbaszadeh和Dehghan使用降阶差分方法对分数阶GS模型进行了数值模拟，但是该方法时间精度仅为一阶。在文 [7] 中，Pindzaa和Owolabi将Fourier谱方法应用于含有分数阶导数的反应–扩散模型进行求解，他们考虑了包括分数阶GS模型和分数阶Schnakenberg系统的几类模型。在文 [8] 中，Alzahrani与Khaliq用Fourier谱方法和高阶时间步进方法模拟了一些空间分数阶反应–扩散模型，包括分数阶Brusselator系统、分数阶Schnakenbergx系统和分数阶GS模型。但是文 [7] [8] 都没有给出全离散数值格式的理论分析。由于分数阶算子具有非局部性，数值上将会获得稠密的或者满的系数矩阵。为了克服计算量过大的瓶颈，Wang提出了一类快速数值计算算法，可以高效地解决这一问题 [9] [10] [11] [12]。樊恩宇利用有限元法研究了GS方程的数值解 [13]。王亭亭采用有限差分方法导出GS方程的数值格式，并对时间半离散数值格式的稳定性进行了分析 [14]。

算子分类方法 [15] [16] 是求解复杂时间依赖模型的有效方法，其主要思想是将原始问题分解为一系列简单的子问题，通过构造子问题的行之有效的格式，达到整体算法的最优，因此算子分裂方法已被广泛应用于解决许多复杂问题。而本文就是基于算子分裂法，提出了求解分数阶GS模型的一种高效数值逼近算法，在时间方向上用算子分裂法将GS方程分解成线性部分和非线性部分，空间方向上采用差分方法。针对线性子问题，本文采用差分方法求解；针对非线性子问题，本文采用Crank-Nicolson差分格式和Rubin-Graves线性化技术处理非线性项。最后通过数值算例验证格式的收敛阶并对长时间动力行为进行模拟。对比前人的研究，本文主要有以下三个创新点：

· 提出了一种求解分数阶GS模型的新型线性化算子分裂格式；

· 分析了格式在L²-norm下的稳定性和收敛性；

· 对比研究了不同 $\alpha$ 对数值实验结果的影响。

2. 一种新的线性化二阶算子分裂方法

2.1. Strang算子分裂方法

我们现在描述一下分裂策略。让 $S_{\mathcal{L}}$ 和 $S_N$ 是与以下线性部分和非线性部分相关的精确解算子。

线性部分为：

$\{\begin{array}{l}{u}_t=-{\mu }_u{\left(-\Delta \right)}^{\alpha /2}u-Fu\\ v_t=-{\mu }_v{\left(-\Delta \right)}^{\alpha /2}v-\left(F+\kappa \right)v\end{array}$ (2.1)

非线性部分为：

$\{\begin{array}{l}{u}_t=-u{v}^2+F\\ v_t=u{v}^2\end{array}$ (2.2)

然后我们利用如下的二阶对称Strang分裂法估计方程(1.1)由t到 $t+\tau$ 的近似解：

$u\left(x,t+\tau \right)\approx S_L\left(\frac{\tau }{2}\right)S_N\left(\tau \right)S_L\left(\frac{\tau }{2}\right)u\left(x,t\right)$ (2.3)

其中 $u={\left(u,v\right)}^{\text{T}}$

接下来，我们将用近似解算子 $S_{\mathcal{L}}^{\tau ,h}$ 和 $S_N^{\tau ,h}$ 近似代替精确解算子 $S_{\mathcal{L}}$ 和 $S_N$ 。

2.2. 线性子问题的数值近似： $S_{\mathcal{L}}\to S_{\mathcal{L}}^{\tau ,h}$

针对线性子问题(2.1)，我们研究整数阶的情况，整数阶情况由如下给出：

$\{\begin{array}{l}{u}_t={\mu }_u{u}_{xx}-Fu\\ v_t={\mu }_v{v}_{xx}-\left(F+\kappa \right)v\end{array}$ (2.4)

接下来我们采用Crank-Nicolson方法对其进行离散化得到：

$\{\begin{array}{l}\frac{u_i^{k+1}-u_i^k}{\tau }=\frac{{\mu }_u}{2}\left(\frac{u_{i+1}^{k+1}-2u_i^{k+1}+u_{i-1}^{k+1}}{h^2}+\frac{u_{i+1}^k-2u_i^k+u_{i-1}^k}{h^2}\right)-\frac{F}{2}\left(u_i^{k+1}+u_i^k\right)\\ \frac{v_i^{k+1}-v_i^k}{\tau }=\frac{{\mu }_v}{2}\left(\frac{v_{i+1}^{k+1}-2v_i^{k+1}+v_{i-1}^{k+1}}{h^2}+\frac{v_{i+1}^k-2v_i^k+v_{i-1}^k}{h^2}\right)-\left(F+\kappa \right)\frac{v_i^{k+1}+v_i^k}{2}\end{array}$ (2.5)

令 $a=\frac{{\mu }_u\tau }{2h^2}$ ， $b=1+\frac{\tau {\mu }_u}{h^2}+\frac{F\tau }{2}$ ， $c=1-\frac{\tau {\mu }_u}{h^2}-\frac{F\tau }{2}$ ， $d=\frac{{\mu }_v\tau }{2h^2}$ ， $e=1+\frac{\tau {\mu }_v}{h^2}+\frac{\left(F+\kappa \right)\tau }{2}$ ，

$f=1-\frac{\tau {\mu }_v}{h^2}-\frac{\left(F+\kappa \right)\tau }{2}$

进一步简化得到：

$\{\begin{array}{l}A{U}^{k+1}=B{U}^k\\ C{V}^{k+1}=D{V}^k\end{array}$ (2.6)

其中：

$A=\left(\begin{array}{ccccc}b& -a& & & \\ -a& b& -a& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & -a\\ & & & -a& b\end{array}\right)$ ， $B=\left(\begin{array}{ccccc}c& a& & & \\ a& c& a& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & a\\ & & & a& c\end{array}\right)$ ， $C=\left(\begin{array}{ccccc}e& -d& & & \\ -d& e& -d& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & -d\\ & & & -d& e\end{array}\right)$

$D=\left(\begin{array}{ccccc}f& d& & & \\ d& f& d& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & d\\ & & & d& f\end{array}\right)$ ， $U^{k+1}=\left(\begin{array}{c}{u}_1^{k+1}\\ ⋮\\ u_{m-1}^{k+1}\end{array}\right)$ ， $U^k=\left(\begin{array}{c}{u}_1^k\\ ⋮\\ u_{m-1}^k\end{array}\right)$ ， $V^{k+1}=\left(\begin{array}{c}{v}_1^{k+1}\\ ⋮\\ v_{m-1}^{k+1}\end{array}\right)$ ， $V^k=\left(\begin{array}{c}{v}_1^k\\ ⋮\\ v_{m-1}^k\end{array}\right)$

为进一步验证格式的有效性，将对其做稳定性分析。

定理1：对于任意 $\tau \ge 0$ ， $h\ge 0$ ，差分格式(2.5)是无条件稳定的。

证明：对于(2.5)式，可以简化为：

$\{\begin{array}{l}-a{u}_{i-1}^{k+1}+b{u}_i^{k+1}-a{u}_{i+1}^{k+1}=a{u}_{i-1}^k+c{u}_i^k+a{u}_{i+1}^k\\ -d{v}_{i-1}^{k+1}+e{v}_i^{k+1}-d{v}_{i+1}^{k+1}=d{v}_{i-1}^k+f{v}_i^k+d{v}_{i+1}^k\end{array}$

利用Fourier方法，令 $u_j^n=E^n{\text{e}}^{i\omega x_j}$ ， $v_j^n=F^n{\text{e}}^{i\omega x_j}$ ，得：

$\{\begin{array}{l}{E}^{n+1}=\frac{1-\frac{F\tau }{2}+a\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}{1+\frac{F\tau }{2}-a\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}{E}^n\\ F^{n+1}=\frac{1-\frac{\left(F+\kappa \right)\tau }{2}+d\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}{1+\frac{\left(F+\kappa \right)\tau }{2}-d\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}{F}^n\end{array}$

由此得到增长因子：

$\begin{array}{c}{G}_u\left(\tau ,\omega \right)=\frac{1-\frac{F\tau }{2}+a\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}{1+\frac{F\tau }{2}-a\left({\text{e}}^{i\omega h}+{\text{e}}^{-i\omega h}-2\right)}\\ =\frac{1-\frac{F\tau }{2}+2a\left(\cos \left(\omega h\right)-1\right)}{1+\frac{F\tau }{2}-2a\left(\cos \left(\omega h\right)-1\right)}\\ =\frac{1-\frac{F\tau }{2}-4a{\sin }^2\frac{\omega h}{2}}{1+\frac{F\tau }{2}+4a{\sin }^2\frac{\omega h}{2}}\end{array}$

从而 $\left|G_u\left(\tau ,\omega \right)\right|\le 1$ ，同理可得： $\left|G_v\left(\tau ,\omega \right)\right|\le 1$ ，从而此格式无条件稳定。

2.3. 非线性子问题的数值近似： $S_N\to S_N^{\tau ,h}$

基于Crank-Nicolson差分格式和Rubin-Graves线性化技术处理 $u{v}^2$ 。周期边界条件下的非线性子问题(2.2)利用以下线性化格式求解：

$\{\begin{array}{l}\frac{u_i^{m+1}-u_i^m}{\tau }=-\frac{1}{2}{\left(v_i^m\right)}^2{u}_i^{m+1}-u_i^m{v}_i^m{v}_i^{m+1}+\frac{1}{2}{\left(v_i^m\right)}^2{u}_i^m+F\\ \frac{v_i^{m+1}-v_i^m}{\tau }=\frac{1}{2}{\left(v_i^m\right)}^2{u}_i^{m+1}+u_i^m{v}_i^m{v}_i^{m+1}-\frac{1}{2}{\left(v_i^m\right)}^2{u}_i^m\end{array}$ (2.6)

现使用“冻结系数”策略来讨论上述格式的稳定性。由以下常数冻结 ${\left(v^m\right)}^2$ 和 $u^m{v}^m$ 两项：

${\theta }_1:=\underset{0\le m\le M}{\max }\underset{0\le i}{\max }\left\{{\left(v_i^m\right)}^2\right\}$ ， ${\theta }_2:=\underset{0\le m\le M}{\max }\underset{0\le i}{\max }\left\{u_i^m{v}_i^m\right\}$

那么(2.6)用矩阵可以表示为：

$X{u}_i^{m+1}=Y{u}_i^m+Z$ (2.7)

其中：

$u_i^m=\left(\begin{array}{c}{u}_i^m\\ v_i^m\end{array}\right)$ ， $X=\left(\begin{array}{cc}1+\frac{\tau }{2}{\theta }_1& \tau {\theta }_2\\ -\frac{\tau }{2}{\theta }_1& 1-\tau {\theta }_2\end{array}\right)$ ， $Y=\left(\begin{array}{cc}1+\frac{\tau }{2}{\theta }_1& 0\\ 0& 1-\frac{\tau }{2}{\theta }_2\end{array}\right)$ ， $Z=\left(\begin{array}{c}\tau F\\ 0\end{array}\right)$

当 $1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\ne 0$ 时，矩阵A可逆。

于是有 $X^{-1}Y=\frac{1}{1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1}\left(\begin{array}{cc}1+\frac{\tau }{2}{\theta }_1-\tau {\theta }_2-\frac{{\tau }^2}{2}{\theta }_1{\theta }_2& -\tau {\theta }_2\left(1-\frac{\tau }{2}{\theta }_2\right)\\ \frac{\tau }{2}{\theta }_1\left(1+\frac{\tau }{2}{\theta }_1\right)& 1-\frac{\tau }{2}{\theta }_2+\frac{\tau }{2}{\theta }_1-\frac{{\tau }^2}{4}{\theta }_1{\theta }_2\end{array}\right)$

记：

$\{\begin{array}{l}\begin{array}{cc}{g}_{11}=1+\frac{\tau }{2}{\theta }_1-\tau {\theta }_2-\frac{{\tau }^2}{2}{\theta }_1{\theta }_2,& g_{12}=-\tau {\theta }_2\left(1-\frac{\tau }{2}{\theta }_2\right)\end{array}\\ \begin{array}{cc}{g}_{21}=\frac{\tau }{2}{\theta }_1\left(1+\frac{\tau }{2}{\theta }_1\right),& g_{22}=1-\frac{\tau }{2}{\theta }_2+\frac{\tau }{2}{\theta }_1-\frac{{\tau }^2}{4}{\theta }_1{\theta }_2\end{array}\end{array}$

于是有 ${\left(X^{-1}Y\right)}^{\text{T}}\left(X^{-1}Y\right)=\frac{1}{{\left(1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\right)}^2}\left(\begin{array}{cc}{g}_{11}^2+g_{21}^2& g_{11}{g}_{12}+g_{21}{g}_{22}\\ g_{11}{g}_{12}+g_{21}{g}_{22}& g_{12}^2+g_{22}^2\end{array}\right)$

假设：

$\tau \le \frac{1}{\left|{\theta }_1-2{\theta }_2\right|}$ (2.8)

于是存在一个常数 $C_1$ ，使得

$\begin{array}{l}\frac{g_{11}^2+g_{21}^2+\left|g_{11}{g}_{12}+g_{21}{g}_{22}\right|}{{\left(1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\right)}^2}\le \frac{g_{11}^2+g_{21}^2+\left|g_{11}{g}_{12}\right|+\left|g_{21}{g}_{22}\right|}{{\left(1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\right)}^2}\\ ={\left(1-\frac{\frac{{\tau }^2}{2}{\theta }_1{\theta }_2}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right)}^2+\frac{{\tau }^2}{4}\frac{{\theta }_1^2{\left(1+\frac{\tau }{2}{\theta }_1\right)}^2}{{\left(1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)\right)}^2}+\tau \left|\frac{{\theta }_2\left(1-\frac{\tau }{2}{\theta }_2\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\cdot \left|1-\frac{\frac{\tau }{2}{\theta }_1{\theta }_2}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\\ +\frac{\tau }{2}\left|\frac{{\theta }_1\left(1+\frac{\tau }{2}{\theta }_1\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\cdot \left|1+\frac{\frac{\tau }{2}{\theta }_2\left(1-\frac{\tau }{2}{\theta }_1\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\\ \le 1+C_1\tau \end{array}$

$\begin{array}{l}\frac{\left|g_{11}{g}_{12}+g_{21}{g}_{22}\right|+g_{12}^2+g_{22}^2}{{\left(1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\right)}^2}\le \frac{\left|g_{11}{g}_{12}\right|+\left|g_{21}{g}_{22}\right|+g_{12}^2+g_{22}^2}{{\left(1-\tau {\theta }_2+\frac{\tau }{2}{\theta }_1\right)}^2}\\ =\tau \left|\frac{{\theta }_2\left(1-\frac{\tau }{2}{\theta }_2\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\cdot \left|1-\frac{\frac{\tau }{2}{\theta }_1{\theta }_2}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|+\frac{\tau }{2}\left|\frac{{\theta }_1\left(1+\frac{\tau }{2}{\theta }_1\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\cdot \left|1+\frac{\frac{\tau }{2}{\theta }_2\left(1-\frac{\tau }{2}{\theta }_1\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right|\\ +{\tau }^2\frac{{\theta }_2^2{(1-\frac{\tau }{2}{\theta }_2)}^2}{{\left(1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)\right)}^2}+{\left(1+\frac{\frac{\tau }{2}{\theta }_2\left(1-\frac{\tau }{2}{\theta }_1\right)}{1+\tau \left(\frac{{\theta }_1}{2}-{\theta }_2\right)}\right)}^2\\ \le 1+C_1\tau \end{array}$

结合上述两个不等式，我们得到

${\Vert X^{-1}Y\Vert }_2^2\le {\Vert {\left(X^{-1}Y\right)}^{\text{T}}\left(X^{-1}Y\right)\Vert }_{\infty }\le 1+C_1\tau$ (2.9)

其中 ${\Vert X\Vert }_2$ 和 ${\Vert X\Vert }_{\infty }$ 是X的光谱范数和各列之和的最大值。

因此，通过(2.7)和(2.9)我们有如下定理：

定理2：在(2.8)的条件下，对于任何一种网格剖分方法 $U=\left\{{\left(u_i,v_i\right)}^{\text{T}}|0\le i\le N-1\right\}$ 有：

$\Vert S_N^{\tau ,h}U\Vert \le \sqrt{1+C_1\tau }\Vert U\Vert$

现在，对于问题(1.1)我们得到一种快速有效的二阶算子分裂格式：

$\{\begin{array}{l}{U}^{*}=A^{-1}B\left(\frac{\tau }{2}\right)U^m,\\ V^{*}=C^{-1}D\left(\frac{\tau }{2}\right)V^m,\\ \frac{U^{**}-U^{*}}{\tau }=-\frac{1}{2}{\left(V^{*}\right)}^2{U}^{**}-U^{*}{V}^{*}{V}^{**}+\frac{1}{2}{\left(V^{*}\right)}^2{U}^{*}+F,\\ \frac{V^{**}-V^{*}}{\tau }=\frac{1}{2}{\left(V^{*}\right)}^2{U}^{**}+U^{*}{V}^{*}{V}^{**}-\frac{1}{2}{\left(V^{*}\right)}^2{U}^{*},\\ U^{m+1}=A^{-1}B\left(\frac{\tau }{2}\right)U^{**},\\ V^{m+1}=C^{-1}D\left(\frac{\tau }{2}\right)V^{**}.\end{array}$ (2.10)

3. 数值实验

由于分数阶GS模型极难求出解析解，本文采用以下公式计算时间误差及收敛阶：

$\{\begin{array}{l}\max u\left(\tau ,h\right)=\max \left|U_h^{\tau }-U_h^{\tau /2}\right|\\ \text{rate\_}u={\log }_2\left(\frac{\max u\left(\tau ,h\right)}{\max u\left(\frac{\tau }{2},h\right)}\right)\end{array},h足够小$

空间误差及收敛阶计算公式如下：

$\{\begin{array}{l}\max u\left(\tau ,h\right)=\max \left|U_h^{\tau }-U_{h/2}^{\tau }\right|\\ \text{rate\_}u={\log }_2\left(\frac{\max u\left(\tau ,h\right)}{\max u\left(\tau ,h/2\right)}\right)\end{array},\tau 足够小$

其中 $U_h^{\tau }$ 代表采用时空步长分别为 $\tau$ 和h，计算出的T时刻的数值解，变量V采用类似方法计算，所有计算均在VivoBook_ASUSLaptop X512JP_V5000JP计算机上使用Matlab2020a编程软件以双倍精度进行。

收敛性测试与一维相图仿真

考虑区域 $-1\le x\le 1,0\le t\le 1$ 上的空间分数阶Gray-Scott方程，初值如下：

$\{\begin{array}{l}u\left(x,0\right)=\left(t^3+1\right)\sin \left(2\pi x\right)\\ v\left(x,0\right)=\left(t^3+1\right){\cos }^5\left(-\pi x\right)\end{array}$ (3.1)

参数选取如下：

$F=0.03$ ， $\kappa =0.062$ ， ${\mu }_u=2\times {10}^{-5}$ ， ${\mu }_v={10}^{-5}$

下面我们计算 $T=1$ 的不同时间步长、不同空间步长时的最大误差 $Er{r}_{\infty }$ 及相应的收敛阶，结果见表1、表2：

Table 1. Maximum error and convergence order for a spatial step of 1/500

表1. 空间步长为1/500的最大误差及收敛阶

Table 2. Maximum error and convergence order for a time step of 1/4000

表2. 时间步长为1/4000的最大误差及收敛阶

由上表可知，随着网格的剖分变细， $Er{r}_{\infty }$ 变得越来越小，且收敛精度越来越接近理论结果。

同时，利用初值条件(3.1)我们进行了一维相图仿真，并模拟了长时间动力行为( $T=100$ )，结果见图1、图2：

Figure 1. Image of u

图1. u的图像

Figure 2. Image of v

图2. v的图像

4. 结论

本文综合利用算子分裂法、Crank-Nicolson差分格式和Rubin-Graves线性化技术求解空间分数阶Gray-Scott方程，建立求解了GS方程的一种高效数值格式，并给出了稳定性分析。最后，通过数值实验验证了所提出格式的收敛阶。

NOTES

^*通讯作者。