=Paper= {{Paper |id=Vol-2748/Paper59 |storemode=property |title=Numerical Solutions of the Variable-Order Space Fractional Activator-Inhibitor Reaction-Diffusion Systems |pdfUrl=https://ceur-ws.org/Vol-2748/IAM2020_paper_59.pdf |volume=Vol-2748 |authors=Redouane Douaifia,Salem Abdelmalek |dblpUrl=https://dblp.org/rec/conf/iam/Douaifia020 }} ==Numerical Solutions of the Variable-Order Space Fractional Activator-Inhibitor Reaction-Diffusion Systems== https://ceur-ws.org/Vol-2748/IAM2020_paper_59.pdf
Numerical solutions of the variable-order space
fractional activator-inhibitor reaction-diffusion
systems
Redouane Douaifiaa , Salem Abdelmaleka,b
a
    Laboratory of Mathematics, Informatics and Systems (LAMIS), Larbi Tebessi University – Tebessa, Algeria
b
    Department of Mathematics and Computer Science, Larbi Tebessi University – Tebessa, Algeria


                                         Abstract
                                         Reaction–diffusion equations containing fractional derivatives can provide adequate mathematical mod-
                                         els for explaining anomalous diffusion and transport dynamics in complex systems that can not be ad-
                                         equately modeled by standard numerical order equations. Researchers have recently found that many
                                         physical processes display fractional order dynamics that differs with time or space. The continuity of
                                         order in the fractional calculation allows the order of the fractional operator to be regarded as a variable.
                                         The Samko–Ross variable-order fractional operator can be viewed as a generalization of the Riemann–
                                         Liouvile type definition. Although this definition is the most appropriate definition having fundamen-
                                         tal characteristics that are desirable for physical modeling, numerical methods for reaction–diffusion
                                         activator–inhibitor systems using this definition have not yet appeared in the literature. In this paper,
                                         we provide a numerical method to get the approximate solutions of the variable-order space-fractional
                                         Activator–inhibitor systems, namely the reaction–diffusion system with cubic nonlinearity and Brus-
                                         selator model on a 1–D space–domain, we adopt the Riesz variable–order space fractional derivative.
                                         Numerical simulations demonstrate that the finite difference approach is computationally efficient.

                                         Keywords
                                         Variable–order calculus, activator–inhibitor systems, numerical simulation




1. Introduction
In recent decades, fractional derivatives have been commonly used in the modeling of complex
physical and mechanical phenomena in wide classes of complex media with hereditary, fractal
and non-markovian properties, is a field of rapidly growing interest with applications in many
different fields. These include the study of biology [1], hydrology [2], biochemistry [3], finance
[4], and physics [5]. But researchers have found that many important dynamic processes show
fractional order behavior that may change with time and / or space. This fact suggests that
differential calculus is a natural candidate to provide an effective mathematical framework to
describe the complex dynamic problems that appear in different biological and chemical mod-
els, such as viscoelastic materials, anomalous diffusion. In reaction diffusion systems which
can be described by two variables, one has the notion of an activator and an inhibitor. In that
case, be the diffusion coefficient of the inhibitor greater than the diffusion coefficient of the
IAM’20: Third conference on informatics and applied mathematics, 21–22 October 2020, Guelma, ALGERIA
" redouane.douaifia@univ–tebessa.dz (R. Douaifia); salem.abdelmalek@univ–tebessa.dz (S. Abdelmalek)
 0000-0001-6960-4572 (R. Douaifia); 0000-0001-9762-9654 (S. Abdelmalek)
                                       © 2020 Copyright for this paper by its authors.
                                       Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
    CEUR
    Workshop
    Proceedings
                  http://ceur-ws.org
                  ISSN 1613-0073       CEUR Workshop Proceedings (CEUR-WS.org)
activator if spontaneous steady-state patterns are to occur. Thus, pattern formation is said to
occur as a result of "short range activation" and "long range inhibition", the activator-inhibitor
reaction-diffusion systems model was used to simulate pattern forming processes. In this work
we show the behavior of solutions for activator-inhibitor models with variable superdiffusion
by remplacing the second derivative in space by the variable order Riesz fractional deriva-
tive. The studies of dynamics and numerical simulations of fractional-space activator–inhibitor
reaction–diffusion systems have recently met with the interest of researchers (cf. [6, 7]).


2. Variable–order fractional derivative
In [8, 9] Samko and Ross directly generalized the Riemann-Liouvile and Marchaud fractional
integration and differentiation of the case of variable order and discussed some properties and
the inversion formula . It has become a controversial research and has raised widespread con-
cern in the last years.
   This section is devoted to the most important definitions of variable–order direvatives.

2.1. Variable–order Riemann–Liouville fractional derivatives
Left derivative:
                                                            𝜉
                𝛼(𝑡,𝑥)                     1        𝜕𝑚
             𝑎 𝐷𝑥      𝑢(𝑡, 𝑥) =                              (𝜉 − 𝜆)𝑚−𝛼(𝑡,𝑥)−1 𝑢(𝑡, 𝜆)𝑑𝜆
                                   [ Γ(𝑚 − 𝛼(𝑡, 𝑥)) 𝜕𝜉 𝑚 ∫𝑎                               ]𝜉 =𝑥

Right derivative:
                                                            𝑏
                𝛼(𝑡,𝑥)                  (−1)𝑚       𝜕𝑚
             𝑥 𝐷𝑏      𝑢(𝑡, 𝑥) =                              (𝜆 − 𝜉 )𝑚−𝛼(𝑡,𝑥)−1 𝑢(𝑡, 𝜆)𝑑𝜆
                                   [ Γ(𝑚 − 𝛼(𝑡, 𝑥)) 𝜕𝜉 𝑚 ∫𝜉                                ]
                                                                                            𝜉 =𝑥

where Γ(.) is gamma function and 𝑚 − 1 < 𝛼(𝑡, 𝑥) ⩽ 𝑚, 𝑚 ∈ ℕ.

2.2. Variable-order Grünwald-Letnikov fractional derivatives
Left derivative:
                                                             𝑛
                   𝐺 𝛼(𝑡,𝑥)                                                𝛼(𝑡, 𝑥)
                   𝑎 𝐷𝑥     𝑢(𝑡, 𝑥) =       lim      ℎ−𝛼(𝑡,𝑥) ∑(−1)𝑗               𝑢(𝑡, 𝑥 − 𝑗ℎ)
                                        ℎ→0,𝑛ℎ=𝑥−𝑎
                                                            𝑗=0
                                                                       (      𝑗 )

Right derivative:
                                                             𝑛
                   𝐺 𝛼(𝑡,𝑥)                                             𝛼(𝑡, 𝑥)
                   𝑥 𝐷𝑏     𝑢(𝑡, 𝑥) =       lim      ℎ−𝛼(𝑡,𝑥) ∑(−1)𝑗            𝑢(𝑡, 𝑥 + 𝑗ℎ)
                                        ℎ→0,𝑛ℎ=𝑏−𝑥
                                                            𝑗=0
                                                                       ( 𝑗 )

  where
                                       𝛼(𝑡, 𝑥)      Γ(𝛼(𝑡, 𝑥) + 1)
                                               =                         ,
                                   (      𝑗 ) Γ(𝑗 + 1)Γ(𝛼(𝑡, 𝑥) − 𝑗 + 1)
and 𝑚 − 1 < 𝛼(𝑡, 𝑥) ⩽ 𝑚, 𝑚 ∈ ℕ.
2.3. Variable–order Riesz fractional derivative
                𝜕 𝛼(𝑡,𝑥) 𝑢(𝑡, 𝑥)               1
                       𝛼(𝑡,𝑥)
                                   =−                   [
                                                               𝛼(𝑡,𝑥)
                                                            𝑎 𝐷𝑥      𝑢(𝑡, 𝑥) +𝑥 𝐷𝑏𝛼(𝑡,𝑥) 𝑢(𝑡, 𝑥)] .
                  𝜕 |𝑥|                 2 cos ( 𝜋𝛼(𝑡,𝑥)
                                                  2 )


3. General Activator-inhibitor model
Consider a 2-species reaction-diffusion system
                          𝜕𝑢                                   𝜕𝑣
                             = 𝐷𝑢 Δ𝑢 + 𝐹 (𝑢, 𝑣),                  = 𝐷𝑣 Δ𝑣 + 𝐺(𝑢, 𝑣)
                          𝜕𝑡                                   𝜕𝑡
Assume there exists a constant steady state 𝐸𝑒𝑞 = (𝑢𝑒𝑞 , 𝑣𝑒𝑞 ),

                                         𝐹 (𝑢𝑒𝑞 , 𝑣𝑒𝑞 ) = 𝐺(𝑢𝑒𝑞 , 𝑣𝑒𝑞 ) = 0

The Jacobian of the previous system is defined by
                                                               𝜕𝐹         𝜕𝐹
                                                      ⎛ 𝜕𝑢 |𝐸𝑒𝑞           𝜕𝑣 |𝐸𝑒𝑞 ⎞
                                            𝐽   𝐽
                                   𝐽 |𝐸𝑒𝑞 = 11 12 = ⎜                            ⎟
                                           (𝐽21 𝐽22 ) ⎜ 𝜕𝐺                𝜕𝐺      ⎟
                                                      ⎝ 𝜕𝑢 |𝐸𝑒𝑞           𝜕𝑣 |𝐸𝑒𝑞 ⎠

   We say that the reaction diffusion system is the inhibitory activator system if the coefficients
of its Jacobian matrix have the following relationships:

                                          𝐽11 𝐽22 < 0,         𝐽12 𝐽21 < 0
  The species that achieve 𝐽𝑖𝑖 > 0, (𝐽𝑖𝑖 < 0) is called activator (inhibitor), respectively (𝑖 = 1, 2).


4. The variable order space fractional Activator-inhibitor
   models
4.1. Variable order space fractional reaction–diffusion model with cubic
     nonlinearity
We consider the reaction-diffusion model with cubic nonlinearity in one-dimensional space, by
replacing the classical spatial differential operators by variable order Riesz fractional analogues
, we obtain the following system

                   𝜕𝑢       𝜕 𝛼(𝑡,𝑥) 𝑢       1
                      = 𝐷𝑢       𝛼(𝑡,𝑥)
                                        + 𝑢 − 𝑢 3 − 𝑣,                 in      [0, 𝑇 ] × [𝑎, 𝑏] ,      (1)
                   𝜕𝑡      𝜕 |𝑥|             3

                  𝜕𝑣       𝜕 𝛼(𝑡,𝑥) 𝑣
                     = 𝐷𝑣             + 𝑢 − 𝑣 + 𝐴,                      in      [0, 𝑇 ] × [𝑎, 𝑏] ,     (2)
                  𝜕𝑡      𝜕 |𝑥|𝛼(𝑡,𝑥)
and the initial and boundary conditions

                           𝑢(0, 𝑥) = 𝑢0 (𝑥), 𝑣(0, 𝑥) = 𝑣0 (𝑥),           on      [𝑎, 𝑏] ,              (3)
                 𝑢(𝑡, 𝑎) = 𝑢(𝑡, 𝑏) = 𝐵𝐶𝑢 ,          𝑣(𝑡, 𝑎) = 𝑣(𝑡, 𝑏) = 𝐵𝐶𝑣 , ∀𝑡 ∈ [0, 𝑇 ] .                     (4)
Here 𝑢, 𝑣 represent the concentrations of two species having the diffusion rates 𝐷𝑢 , 𝐷𝑣 > 0,
𝐴 ∈ ℝ, 𝐵𝐶𝑢 ⩾ 0 and 𝐵𝐶𝑣 ⩾ 0 are external parameters.

4.2. Variable order space fractional Brusselator model
We consider the well-known reaction-diffusion model Brusselator in one-dimensional space, by
replacing the classical spatial differential operators by variable order Riesz fractional analogues
, we get the following system

               𝜕𝑢       𝜕 𝛼(𝑡,𝑥) 𝑢
                  = 𝐷𝑢             + 𝐴 − (𝐵 + 1)𝑢 + 𝑣𝑢 2 ,                  in        [0, 𝑇 ] × [𝑎, 𝑏] ,         (5)
               𝜕𝑡      𝜕 |𝑥|𝛼(𝑡,𝑥)
             𝜕𝑣       𝜕 𝛼(𝑡,𝑥) 𝑣
                = 𝐷𝑣             + 𝐵𝑢 − 𝑣𝑢 2 ,                                 in           [0, 𝑇 ] × [𝑎, 𝑏] ,   (6)
             𝜕𝑡      𝜕 |𝑥|𝛼(𝑡,𝑥)
and the initial and boundary conditions
                         𝑢(0, 𝑥) = 𝑢0 (𝑥), 𝑣(0, 𝑥) = 𝑣0 (𝑥),            on          [𝑎, 𝑏] ,                     (7)
                 𝑢(𝑡, 𝑎) = 𝑢(𝑡, 𝑏) = 𝐵𝐶𝑢 ,          𝑣(𝑡, 𝑎) = 𝑣(𝑡, 𝑏) = 𝐵𝐶𝑣 , ∀𝑡 ∈ [0, 𝑇 ] .                     (8)
  The species 𝑢, 𝑣 represent the concentrations of two intermediary reactants having the dif-
fusion rates 𝐷𝑢 , 𝐷𝑣 > 0 and 𝐴, 𝐵 > 0 are fixed concentrations, and 𝐵𝐶𝑢 , 𝐵𝐶𝑣 are external
parameters (nonnegative numbers).


5. Explicit Euler Approximation
We show the approximate solutions for systems (1)-(4) and (5)-(8) by applying the explicit finite
difference method which described in [10].
   We consider the numerical approximation in the time domain [0, 𝑇 ] and the space domain
[𝑎, 𝑏]. Let 𝑡𝑘 = 𝑘Δ𝑡 (0 ⩽ 𝑡𝑘 ⩽ 𝑇 ), 𝑘 = 0, … , 𝑀,
   𝑥𝑖 = 𝑎 + 𝑖Δ𝑥 (𝑎 ⩽ 𝑥𝑖 ⩽ 𝑏), 𝑖 = 0, … , 𝑁 , where the time step is Δ𝑡 = 𝑀
                                                                         𝑇
                                                                           and the space step is
                                                                               sec(𝜋𝛼 𝑘 )
      𝑁 . We denote that 𝑢𝑖 = 𝑢(𝑡𝑘 , 𝑥𝑖 ), 𝛼𝑖 = 𝛼(𝑡𝑘 , 𝑥𝑖 ), 𝑐𝑖 = −
Δ𝑥 = 𝑏−𝑎                                                                . The approximate Grün-
                           𝑘                𝑘                 𝑘       𝑖
                                                                    2
wald formulas for the variable-order Riemann-Liouville fractional derivatives approximation
show as follow:
                                                                             𝑖+1
                         𝛼𝑘𝑖             𝐺    𝑖𝛼𝑘                 −𝛼𝑖+1 𝑘  (𝑗) 𝑘
                      𝑎 𝐷𝑥𝑖 𝑢(𝑡𝑘 , 𝑥𝑖 ) =𝑎 𝐷𝑥𝑖 𝑢(𝑡𝑘 , 𝑥𝑖 ) ≈ (Δ𝑥)       ∑ 𝑔𝛼 𝑘 𝑢𝑖+1−𝑗
                                                                                      𝑖+1
                                                                             𝑗=0
                                                                            𝑁 −𝑖+1
                        𝛼𝑖𝑘            𝐺   𝛼𝑖𝑘                    𝑘
                                                                −𝛼𝑖−1                𝑘(𝑗)
                    𝑥𝑖 𝐷𝑏 𝑢(𝑡𝑘 , 𝑥𝑖 ) =𝑥𝑖 𝐷𝑏 𝑢(𝑡𝑘 , 𝑥𝑖 ) ≈ (Δ𝑥)              ∑ 𝑔𝛼 𝑘 𝑢𝑖−1+𝑗
                                                                                       𝑖−1
                                                                             𝑗=0

where 𝑔𝛼(𝑗)𝑘 is the Grünwald weights defined by
         𝑖

                                                        𝛼𝑖𝑘 + 1
                        𝑔𝛼(0)𝑘 = 1, 𝑔𝛼(𝑗)𝑘 =       1−           𝑔 (𝑗−1) , (𝑗 = 1, 2, … )
                          𝑖            𝑖       (           𝑗 ) 𝛼𝑖𝑘
Therefore, the equations (1)-(2) and (5)-(6) can be discretized as follow:
                                      𝑖+1                    𝑁 −𝑖+1
            𝑢𝑖𝑘+1 = 𝑢𝑖𝑘 + 𝐷𝑢     (1)
                                𝑟𝑖,𝑘     (𝑗)   𝑘
                                     ∑ 𝑔𝑖+1,𝑘 𝑢𝑖+1−𝑗
                                                        (2)
                                                     + 𝑟𝑖,𝑘     (𝑗)   𝑘
                                                            ∑ 𝑔𝑖−1,𝑘 𝑢𝑖−1+𝑗   + Δ𝑡𝐹 (𝑢𝑖𝑘 , 𝑣𝑖𝑘 )           (9)
                               (     𝑗=0                    𝑗=0             )
                                      𝑖+1                    𝑁 −𝑖+1
           𝑣𝑖𝑘+1 = 𝑣𝑖𝑘 + 𝐷𝑣     𝑟 (1) ∑ 𝑔 (𝑗) 𝑣𝑖+1−𝑗
                                               𝑘        (2)
                                                     + 𝑟𝑖,𝑘     (𝑗)   𝑘
                                                            ∑ 𝑔𝑖−1,𝑘 𝑣𝑖−1+𝑗   + Δ𝑡𝐺(𝑢𝑖𝑘 , 𝑣𝑖𝑘 )           (10)
                               ( 𝑖,𝑘 𝑗=0 𝑖+1,𝑘              𝑗=0             )

  where 𝑟𝑖,𝑘 = Δ𝑡𝑐𝑖𝑘 (Δ𝑥)−𝛼𝑖+1 , 𝑟𝑖,𝑘 = Δ𝑡𝑐𝑖𝑘 (Δ𝑥)−𝛼𝑖−1 , 𝑔𝑖,𝑘 = 𝑔𝛼(𝑗)𝑘 and
         (1)                   𝑘  (2)                   𝑘  (𝑗)
                                                                      𝑖


       ⎧
       ⎪ 𝐹 (𝑢 𝑘 , 𝑣 𝑘 ) = 𝑢𝑖𝑘 − 13 (𝑢𝑖𝑘 )3 − 𝑣𝑖𝑘        ⎧
                                                        ⎪ 𝐹 (𝑢 𝑘 , 𝑣 𝑘 ) = 𝐴 − (𝐵 + 1)𝑢𝑖𝑘 + 𝑣𝑖𝑘 (𝑢𝑖𝑘 )2
       ⎪ 𝑖 𝑖                                            ⎪ 𝑖 𝑖
       ⎨                                           𝑜𝑟   ⎨
       ⎪      𝑘 𝑘                                       ⎪
       ⎪
       ⎩𝐺(𝑢𝑖 , 𝑣𝑖 ) =        𝑢𝑖𝑘 − 𝑣𝑖𝑘 + 𝐴              ⎪      𝑘 𝑘
                                                        ⎩𝐺(𝑢𝑖 , 𝑣𝑖 ) =          𝐵𝑢𝑖𝑘 − 𝑣𝑖𝑘 (𝑢𝑖𝑘 )2

  The stability and convergence of the explicit Euler approximation are descused in [10].


6. Numerical experiments
In this section, we show approximate solutions for the systems (1)-(4) and (5)-(8) to demon-
strate the changes in solutions behaviour that arise when the exponent is varied from integer
order to fractional order to variable order, and to identify the differences between solutions.
The computer algorithm for numerical method (9)-(10), was written in Matlab, throughout the
simulations we took the following values: 𝑎 = 0, 𝑏 = 4, 𝑇 = 20, spatial and time steps respec-
                            2
tively Δ𝑥 = 204
                , Δ𝑡 = (Δ𝑥)
                         2 − 0.001.


6.1. Reaction-diffusion model with cubic nonlinearity
We take, 𝐴 = −0.1, 𝐷𝑢 = 0.05, 𝐷𝑣 = 1, 𝐵𝐶𝑢 = 0.0503, 𝐵𝐶𝑣 = 0.07504 and the initial conditions:
                    {
                      𝑢0 (𝑥) = 0.0503 − 10−3 cos(2𝜋𝑥) + 10−3 cos(5𝜋𝑥)
                      𝑣0 (𝑥) = 0.02504 + 10−3 cos(2𝜋𝑥) + 10−3 cos(5𝜋𝑥)

The figure 1 shows the approximate solution for system (1)-(4) for 𝛼(𝑡, 𝑥) = 2, when the deriva-
tive is fractional order 𝛼(𝑡, 𝑥) = 1.35 we can see the sulotions in figure 2. The behaviour of the
solution is particularly interesting for the case 𝛼(𝑡, 𝑥) = 1.5 + 0.4 cos(𝑡𝑥) ∗ sin(2𝑡), see figure 3.
Figure 1: Approximate solution of (1)-(4), for 𝛼(𝑡, 𝑥) = 2, where 𝑢1 = 𝑢 and 𝑢2 = 𝑣.




Figure 2: Approximate solution of (1)-(4), for 𝛼(𝑡, 𝑥) = 1.35, where 𝑢1 = 𝑢 and 𝑢2 = 𝑣.


6.2. Brusselator model
We take, 𝐴 = 1, 𝐵 = 3, 𝐷𝑢 = 0.05, 𝐷𝑣 = 1, 𝐵𝐶𝑢 = 1.09, 𝐵𝐶𝑣 = 2.02504 and the initial conditions:
                     {
                       𝑢0 (𝑥) =     1.09 − 10−2 cos(2𝜋𝑥) + 10−3 cos(5𝜋𝑥)
                       𝑣0 (𝑥) = 2.02504 + 10−2 cos(2𝜋𝑥) + 10−3 cos(5𝜋𝑥)
 Figure 4 shows the behavior of the numerical solution for system (5)-(8) with 𝛼(𝑡, 𝑥) = 2,
when the derivative is fractional order 𝛼(𝑡, 𝑥) = 1.35 we can see the sulotions in figure 5. The
Figure 3: Approximate solution of (1)-(4), for 𝛼(𝑡, 𝑥) = 1.5 + 0.4 cos(𝑡𝑥) ∗ sin(2𝑡), where 𝑢1 = 𝑢 and
𝑢2 = 𝑣.




Figure 4: Approximate solution of (5)-(8), for 𝛼(𝑡, 𝑥) = 2, where 𝑢1 = 𝑢 and 𝑢2 = 𝑣.


behaviour of the solution is particularly interesting for the case 𝛼(𝑡, 𝑥) = 1.5 + 0.4 cos(𝑡𝑥) ∗
sin(2𝑡), see figure 6.
Figure 5: Approximate solution of (5)-(8), for 𝛼(𝑡, 𝑥) = 1.35, where 𝑢1 = 𝑢 and 𝑢2 = 𝑣.




Figure 6: Approximate solution of (5)-(8), for 𝛼(𝑡, 𝑥) = 1.5 + 0.4 cos(𝑡𝑥) sin(2𝑡), where 𝑢1 = 𝑢 and 𝑢2 = 𝑣.


7. Conclusions
In this work, we have got interesting behavior of solutions for activator-inhibitor models with
variable superdiffusion by replacing the classical second derivative in space by the variable
order Riesz fractional derivative of order 1 < 𝛼(𝑡, 𝑥) ⩽ 2, it seems that our numerical results will
open horizons for analytical studies about such types of models and guide them by guessing.
References
 [1] S. B. Yuste, K. Lindenberg, Subdiffusion–limited reactions, Chemical physics 284 (2002)
     169–180. doi:10.1016/S0301-0104(02)00546-3.
 [2] D. A. Benson, S. W. Wheatcraft, M. M. Meerschaert, Application of a fractional advection–
     dispersion equation, Water resources research 36 (2000) 1403–1412. doi:10.1029/
     2000WR900031.
 [3] S. B. Yuste, K. Lindenberg, Reaction front in an a+ b→ c reaction–subdiffusion process,
     Physical Review E 69 (2004) 036126. doi:10.1103/PhysRevE.69.036126.
 [4] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial
     data: an empirical study, Physica A: Statistical Mechanics and its Applications 314 (2002)
     749–755. doi:10.1016/S0378-4371(02)01048-8.
 [5] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional
     dynamics approach, Physics reports 339 (2000) 1–77. doi:10.1016/S0370-1573(00)
     00070-3.
 [6] Z. Hammouch, T. Mekkaoui, F. Belgacem, Numerical simulations for a variable order
     fractional schnakenberg model, in: AIP Conference Proceedings, American Institute of
     Physics, 2014, pp. 1450–1455. doi:10.1063/1.4907312.
 [7] Z. Prytula, Mathematical modelling of nonlinear dynamics in activator–inhibitor systems
     with superdiffusion, Bulletin of the National University of Lviv Polytechnic. Computer
     science and information technology 826 (2015) 230–237.
 [8] S. G. Samko, B. Ross, Integration and differentiation to a variable fractional order, Integral
     transforms and special functions 1 (1993) 277–300. doi:10.1080/10652469308819027.
 [9] S. G. Samko, Fractional integration and differentiation of variable order, Analysis Math-
     ematica 21 (1995) 213–236. doi:10.1007/bf01911126.
[10] P. Zhuang, F. Liu, V. Anh, , I. Turner, Numerical methods for the variable-order fractional
     advection–diffusion equation with a nonlinear source term, SIAM Journal on Numerical
     Analysis 47 (2009) 1760–1781. doi:10.1137/080730597.