A parallelized computational model for multidimensional systems of coupled nonlinear fractional hyperbolic equations

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.1 (1-19)

Journal of Computational Physics ••• (••••) ••••••

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A parallelized computational model for multidimensional systems of coupled nonlinear fractional hyperbolic equations J.E. Macías-Díaz Departamento de Matemáticas y Física, Universidad Autónoma de Aguascalientes, Avenida Universidad 940, Ciudad Universitaria, Aguascalientes, Ags. 20131, Mexico

a r t i c l e

i n f o

Article history: Received 16 July 2019 Received in revised form 11 October 2019 Accepted 14 October 2019 Available online xxxx Keywords: Generalized systems of hyperbolic fractional equations Riesz Space-fractional derivatives Discrete fractional energy method Stability and convergence analysis Parallel computational implementation Complex pattern formation

a b s t r a c t In this work, we consider a general multidimensional system of hyperbolic partial differential equations with fractional diffusion of the Riesz type, constant damping and coupled nonlinear reaction terms. The system generalizes many particular models from the physical sciences (including inhibitor-activator models in chemistry, diffusive nonlinear systems in population dynamics and relativistic wave equations), and considers the presence of an arbitrary number of both spatial dimensions and dependent variables. Motivated by the wide range of applications, we propose an explicit four-step finitedifference methodology to approximate the solutions of the continuous system. The properties of stability, boundedness and convergence of the scheme are proved rigorously using a discrete form of the fractional energy method. An efficient computational implementation of the scheme is also proposed in this work. It is important to recall that algorithms for space-fractional systems are computationally highly demanding. To alleviate this problem, a parallel implementation of our scheme is proposed using a vector reformulation of the numerical method. We provide some illustrative simulations on the formation of complex patterns in the two-dimensional scenario, and even in the computationally intense three-dimensional case. For the sake of convenience, an algorithmic presentation of our computational model is provided in this manuscript. © 2019 Elsevier Inc. All rights reserved.

1. Introduction The investigation on the conditions under which Turing patterns appear in physical systems has been a highly transited avenue of research from the mathematical, the numerical and the physical points of view. It is well known that many nonlinear systems exhibit Turing patterns under suitable conditions on the parameters of the model and the initial data. Some diffusion-reaction systems exhibit the presence of patterns, like some coupled systems which describe the interaction between inhibitor and activator substances in chemistry [1,2], especially in chemical reactions with chlorine dioxide, iodine, and malonic acid [3]. In fact, several different Turing patterns are also found in the chlorine-iodide-malonic acid reaction, including hexagonal, striped [4], and oscillatory structures [5]. Outside the chemical sciences, there is also evidence of the presence of this kind of nonlinear behavior. For instance, there is experimental evidence on the presence of azimuthal Turing patterns in Kerr combs generated by whispering-gallery-mode resonators [6], and the theory suggests that diffusion-reaction systems may be used to understand the formation of this type of patterns in biology [7]. Even the labyrinthine structure

E-mail address: [email protected]. https://doi.org/10.1016/j.jcp.2019.109043 0021-9991/© 2019 Elsevier Inc. All rights reserved.

JID:YJCPH AID:109043 /FLA

2

[m3G; v1.261; Prn:4/11/2019; 13:14] P.2 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

of the cerebral cortex has been identified as a complex three-dimensional Turing pattern [8], among other applications [9]. Moreover, it is well known that the Brusselator is a system which exhibits the presence of that type of structures [10,11]. The physical models described in the paragraph above are diffusive equations. It is well known that these systems have the physical limitation that any perturbation is instantaneously propagated throughout the system. For that reason, hyperbolic forms of these equations are physically preferred. It turns out that a wide variety of Turing patterns have been found also in hyperbolic systems. For example, Turing instabilities have been investigated in diffusive predator-prey models with hyperbolic mortality [12,13], in hyperbolic models for locally interacting cell systems [14], in the hyperbolic chaos of standing-wave patterns generated by a modulated pump source [15] and in the spreading of infectious diseases in hyperbolic susceptible-infected-removed models [16]. Also, the presence of various wave and Turing patterns has been studied in the context of hyperbolic forms of the Brusselator [17,18], in some hyperbolic models for self-organized biological aggregations and movement [19], in network systems through collective patterns and single differentiated nodes [20], in predator-prey reaction-diffusion models with spatio-temporal delays [21] and in hyperbolic vegetation models for semiarid environments [22], just to mention some systems. In recent years, fractional derivatives have been introduced to mathematical models in order to provide more realistic descriptions of the physical phenomena. For instance, many fractional systems have been obtained as the continuous limit of discrete systems of particles with long-range interactions [23,24]. However, independently of that, fractional derivatives have been successfully used in the theory of viscoelasticity [25], the theory of thermoelasticity [26], financial problems under a continuous time frame [27], self-similar protein dynamics [28] and quantum mechanics [29]. Moreover, some distributedorder fractional diffusion-wave equations are used in the modeling of groundwater flow to and from wells [30,31]. As expected, the complexity of fractional problems is considerably higher than that of integer-order models, whence the need to design reliable numerical techniques to approximate the solutions is pragmatically justified [32]. In this direction, the literature reports on various methods to approximate the solutions of fractional systems. For example, some numerical methods have been proposed to solve fractional partial differential equations using fractional centered differences [33,34], the time-fractional diffusion equation [35], the fractional Schrödinger equation in multiple spatial dimensions [36], the nonlinear fractional Korteweg–de Vries–Burgers equation [37], the fractional FitzHugh–Nagumo monodomain model in two spatial dimensions [38], distributed-order time-fractional diffusion-wave equations on bounded domains [39] and some Hamiltonian hyperbolic fractional differential equations that generalize various well-known equations from quantum field theory [40]. In light of these facts, the investigation on the presence of Turing patterns in fractional systems (in both the parabolic and hyperbolic types) has been also an interesting topic of research in recent years. Some reports have studied the presence of these patterns in simple models with fractional diffusion [41], in fractional reaction-diffusion systems with indices of different order [42] and with multiple homogeneous states [43], in hyperbolic inhibitor-activator systems with fractional diffusion [44], in two-species fractional reaction-diffusion systems [45] and in reaction models with sub-diffusion [46]. At the same time, various numerical methods have been proposed in the literature to investigate complex nonlinear models [47–49]. However, it is worth pointing out that most of the computational methods proposed in the literature are computationally slow and highly demanding [50]. This is a consequence of the fact that the calculation of fractional derivatives require the use of global information of the system at each time step [51–54]. The present work is intended to alleviate partially such shortcomings. Moreover, it is important to point out that, as opposed to some previous efforts by this author [33,40,55], the present approach is not Hamiltonian and considers a family consisting of more than one hyperbolic partial differential equation. As a consequence, various physical systems can be simulated using the numerical technique reported in this work. The purpose of this work is to propose a parallel (and, thus, computationally economic [56]) methodology to investigate general coupled hyperbolic systems with Riesz space-fractional diffusion. In order to generalize our physical model as much as possible, the system investigated in this work will consider an arbitrary finite number p of spatial dimensions, as well as an arbitrary finite number q of dependent variables. We will employ fractional centered differences to approximate the fractional partial derivatives of the model in view of their mathematical and computational simplicity, and an explicit fourstep finite-difference discretization of our physical model will be proposed. Obviously, the explicit nature of the scheme will is an advantage when considering a high number of spatial variables. We will establish rigorously the main numerical features of our scheme, namely, its second-order consistency, stability, boundedness and quadratic convergence. The computational implementation of the scheme will be carried out using a convenient reformulation of the numerical method using vector algebra when approximating the fractional derivatives. Various computer languages have efficient implementations of vector algebra operations, including Fortran, C++ and Matlab. Our implementation is carried out using Gfortran on a Linux operating system, and we will increase the speed of the calculations using the multi-processing package OpenMP. As an illustrative application of our computer code, we will exhibit the presence of Turing patterns in some two- and threedimensional inhibitor-activator systems. However, we will note that the methodology proposed in this work can be applied to a wide variety of physical models. This work is organized as follows. In Section 2, we provide the definition of Riesz partial fractional partial derivative, and introduce the mathematical system of interest. We note therein that various models from physics are particular cases of our model, and we close that stage providing a vector form of the mathematical model, considering some parameter simplifications. In Section 3, we introduce the discrete nomenclature and operators which will be used throughout this manuscript, including the concept of fractional centered differences, which is crucial in our investigation. An explicit four-

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.3 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

3

step finite-difference scheme to solve the continuous problem of interest is introduced in that stage. The discrete model is then analyzed rigorously in Section 4. Concretely, we establish that the scheme is consistent of the second order, and we use a discrete form of the fractional energy method to show that the model is stable and quadratically convergent under appropriate parameter conditions. Section 6 shows some illustrative simulations in the two- and the three-dimensional scenarios. Concretely, we exhibit the formation of Turing patterns in those systems. Finally, we close this work with a section of concluding remarks. 2. Preliminaries Throughout this work, we will use the symbol In to represent the set of indexes {1, . . . , n}, and define I n = In ∪ {0}, for each n ∈ N . Moreover, if x ∈ R p then we convey that x = (x1 , . . . , x p ), for any p ∈ N . The present section is devoted to introduce the continuous nomenclature and the mathematical model under investigation in this manuscript. In particular, the following concepts are standard in the literature. They are the cornerstone to describe our continuous model, and they are recalled here for the sake of convenience. Throughout,  denotes the usual Gamma function. Definition 2.1 (Podlubny [57]). Let f : R → R be a function, and let n ∈ N ∪ {0} and fractional derivative of f of order α at x ∈ R is defined (when it exists) as

−1 dn = π α α d|x| 2 cos( 2 )(n − α ) dxn

dα f (x)

∞ −∞

α ∈ R satisfy n − 1 < α < n. The Riesz

f (ξ )dξ . |x − ξ |α +1−n

(2.1)

Definition 2.2 (Podlubny [57]). Assume that p ∈ N and i ∈ I p . Let α ∈ R and n ∈ N satisfy n − 1 < α < n, and suppose that u : R p × R → R is a function. Then the left and the right Riemann–Liouville fractional partial derivative of u of order α with respect to xi at the point (x, t ) ∈ R p × R are defined, respectively, by

−∞

D α u (x, t ) = xi

α

xi D ∞ u (x, t )

=

∂n (n − α ) ∂ xni 1

n

n

(−1) ∂ (n − α ) ∂ xni

x

(xi − ξ )n−1−α u (x1 , . . . , xi −1 , ξ, xi +1 , . . . , x p , t )dξ,

(2.2)

−∞ ∞

(ξ − xi )n−1−α u (x1 , . . . , xi −1 , ξ, xi +1 , . . . , x p , t )dξ,

(2.3)

x

whenever they exist. For the remainder of this work, we will consider differentiation orders satisfying 1 < α < 2. In that case, we define the Riesz partial fractional derivative of u of order α with respect to xi at the point (x, t ) as

 α  ∂ α u (x, t ) 1 α =− a D xi u (x, t ) + xi D b u (x, t ) . α ∂|xi | 2 cos(πα /2)

(2.4)

Moreover, for convenience, we will agree that the Riesz partial fractional derivative of u of order 1 and 2 with respect to xi coincide with the usual first- and second-order partial derivative of u with respect to xi , respectively. Finally, the Riesz fractional Laplacian of u of order α at (x, t ) will be defined as

α u (x, t ) =

p  ∂ α u (x, t ) i =1

∂|xi |α

(2.5)

.

For the remainder of this paper, we let p ∈ N and assume that ai , b i ∈ R satisfy ai < b i , for each i ∈ I p . Let T ∈ R+ , and fix the spatial and the spatial-temporal domains

B=

p  (ai , bi ),

 = B × (0, T ),

(2.6)

i =1

respectively. Obviously, the sets B and  are open in R p +1 , their respective closures will be denoted by B and , and we will use the notation ∂ B to represent the boundary of B. Also, we will fix q ∈ N , and we will suppose that u j :  → R is a function, for each j ∈ Iq . For convenience, we will extend the definition of u j to all the space R p × [0, T ], by letting u j (x, t ) = 0 for each (x, t ) ∈ (R p \ B ) × [0, T ] and j ∈ Iq . Let F j : Rq → R be a function for each j ∈ Iq , and define u :  → Rq by

u (x, t ) = (u 1 (x, t ), u 2 (x, t ), . . . , u q (x, t )),

∀(x, t ) ∈ .

(2.7)

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.4 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

4

Moreover, we will agree that γ j is a nonnegative constant, τ j and d j are positive, α j ∈ (1, 2] and φ j , ψ j : B → R are smooth functions, for each j ∈ Iq . We will investigate computationally the solution of the following problem governed by a system of hyperbolic space-fractional partial differential equations with nonlinear reaction terms:

τj

∂ 2 u j (x, t ) ∂ u j (x, t ) + γj = d j α j u j (x, t ) − F j (u (x, t )), ∀(x, t ) ∈ , ∀ j ∈ Iq , ∂t2 ⎧∂ t u ( x , 0 ) = φ ( x ), ∀ x ∈ B , ∀ j ∈ I , q j ⎪ ⎨ j ∂ u j (x, 0) such that = ψ j (x), ∀x ∈ B , ∀ j ∈ Iq , ⎪ ∂t ⎩ u j (x, t ) = 0, ∀(x, t ) ∈ ∂ B × [0, T ], ∀ j ∈ Iq .

(2.8)

More precisely, the system of partial differential equations in (2.8) can be explicitly described as

⎧ ∂ 2 u 1 (x, t ) ∂ u 1 (x, t ) ⎪ ⎪ τ1 + γ1 = d1 α1 u (x, t ) − F 1 (u (x, t )), ⎪ 2 ⎪ ∂ t ∂ t ⎪ ⎪ 2 ⎪ ⎪ ⎨ τ ∂ u 2 (x, t ) + γ ∂ u 2 (x, t ) = d α2 u (x, t ) − F (u (x, t )), 2 2 2 2 2 ∂t ∂t2 ⎪ . ⎪ .. ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎪ ⎩ τ ∂ uq (x, t ) + γ ∂ uq (x, t ) = d αq u (x, t ) − F (u (x, t )), q q q q q ∂t ∂t2

∀(x, t ) ∈ , ∀(x, t ) ∈ ,

(2.9)

∀(x, t ) ∈ .

Example 2.3. It is important to point out that (2.8) generalizes various mathematical models used in the physical sciences. The following are some examples where (2.8) finds interesting applications. (a) Suppose that q = 1 in model (2.8), and let u = u 1 . If F 1 (u (x, t )) = sin(u (x, t )) then the resulting model is a fractional form of the damped sine-Gordon equation [58,59], which is well-known in mathematical physics. (b) The system (2.8) is a generalization of a predator-prey system with Michaelis–Menten-type functional response [60]. In that model, q = 2, and u 1 (x, t ) and u 2 (x, t ) represent the prey and the predator densities, respectively. After a scaling analysis, that system is described by (2.8) with reactions



F 1 (u 1 , u 2 ) = Ru 1 1 − F 2 (u 1 , u 2 ) =

Su 1 u 2 u 2 + Su 1

u1 S

−

Su 1 u 2 u 2 + Su 1

− Q u2 ,

,

∀ u 1 , u 2 ∈ R+ ,

(2.10)

∀ u 1 , u 2 ∈ R+ ,

(2.11)

where the constants R, S and Q are positive. (c) The model (2.8) is also a hyperbolic and fractional generalization of some systems which describe the interaction between some activator and inhibitor substances in chemistry [1]. Such is the case when q = 2 and

F 1 (u 1 , u 2 ) = u 1 − au 2 + bu 1 u 2 − u 31 ,

(2.12)

F 2 (u 1 , u 2 ) = u 1 − cu 2 .

(2.13)

In that context, u 1 and u 2 represent the amounts of the activator and inhibitor substances, respectively.

2

Like these examples, there are various other models in the physical sciences which are generalized by (2.8) in various ways and, for that reason, its numerical investigation is an important topic of research. For simplification purposes and for the remainder of this work, convey that τ j = 1, γ = γ j ∈ R+ ∪ {0}, d = d j ∈ R+ and α = α j ∈ (1, 2], for each j ∈ Iq . It is worth pointing out that our analysis would be similar for the general scenario, but we have chosen these simplifications for the sake of easiness. Moreover, we will observe the vector notation

 ∂ 2 uq (x, t ) ∂ 2 u 1 (x, t ) ∂ 2 u 2 (x, t ) , ,..., , ∀(x, t ) ∈ , ∂t2 ∂t2 ∂t2   ∂ uq (x, t ) ∂ u 1 (x, t ) ∂ u 2 (x, t ) ∂ u (x, t ) , ∀(x, t ) ∈ , = , ,..., ∂t ∂t ∂t ∂t   α u (x, t ) = α u 1 (x, t ), α u 2 (x, t ), . . . , α uq (x, t ) , ∀(x, t ) ∈ ,   ∀(x, t ) ∈ . F (x, t ) = F 1 (u (x, t )), F 2 (u (x, t )), . . . , F q (x, t ) , ∂ 2 u (x, t ) = ∂t2

(2.14) (2.15) (2.16) (2.17)

Also, let φ(x) = (φ1 (x), φ2 (x), . . . , φq (x)) and ψ(x) = (ψ1 (x), ψ2 (x), . . . , ψq (x)), for each x ∈ B. With this notation and the simplifications proposed above, the mathematical model (2.8) can be expressed alternatively in vector form as

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.5 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

∂ 2 u (x, t ) ∂ u (x, t ) +γ = dα u (x, t ) − F (u (x, t )), ∀(x, t ) ∈ , 2 ∂t ⎧∂ t u (x, 0) = φ(x), ∀x ∈ B , ⎪ ⎨ ∂ u (x, 0) such that = ψ(x), ∀x ∈ B , ⎪ ∂t ⎩ u (x, t ) = 0, ∀(x, t ) ∈ ∂ B × [0, T ].

5

(2.18)

Finally, we will observe the following assumption for the remainder of this work: (H) There exists a smooth function G = (G 1 , . . . , G q ) : Rq → Rq with G 1 , . . . , G q : Rq → R, such that

F j (u ) =

∂ G j (u ) , ∂u j

∀ j ∈ Iq , ∀u ∈ Rq .

(2.19)

As we will note in Section 6, the models described in Example 2.3 satisfy this hypothesis. Moreover, it is also important to mention that a system (2.18) which satisfies this assumption need not be Hamiltonian. This fact implies that the type of systems satisfying (H) is substantially large. As expected, the finite-difference scheme proposed in the following section will be also general enough to consider not necessarily Hamiltonian models, like those nonlinear systems in which complex patterns appear. 3. Numerical model The present section will be devoted to describe the numerical method to approximate the solutions of (2.18). Agree that h i and τ are positive step-sizes, for each i ∈ I p . Let N = T /τ and M i = (b i − ai )/h i be natural numbers, for each i ∈ I p . Consider uniform partitions of [ai , b i ] and [0, T ] given by

xi ,ki = ai + ki h i , t n = nτ ,

∀i ∈ I p , ∀ki ∈ I M i ,

(3.1)

∀n ∈ I N ,

t n + 1 = t n + τ /2 ,

(3.2)

∀n ∈ I N −1 .

2

(3.3)

We introduce the discrete sets

K=

p 

K=

I M i −1 ,

i =1

p 

I Mi ,

(3.4)

i =1

and define xk = (x1,k1 , . . . , x p ,k p ) for each multi-index k = (k1 , . . . , k p ) ∈ K. Let ∂ K represent the boundary of the mesh K, that is, let ∂ K be the set of multi-indexes k ∈ K such that xk ∈ ∂ B. In this manuscript, we will convey that

h = (h1 , . . . , h p ),

(3.5)

h∗ = h1 · · · h p ,

(3.6)

and Rh will represent the spatial mesh {xk : k ∈ K }. Moreover, Vh will denote the vector space of all real functions define on Rh , which vanish at those points on the boundary of B. For simplification purposes, we will agree that w k = w (xk ), for each w ∈ Vh and k ∈ K. For each j ∈ Iq and (k, n) ∈ K × I N , let unj,k = u j (xk , tn ), and let unj,k represent an approximation to the exact value of unj,k . If w ∈ Vh then we will use the notation w = (wk )k∈K . As a consequence, note that unj ∈ Vh can be represented as (unj,k )k∈K , for each j ∈ Iq . Moreover, we will use the convention u j = (unj )n∈I

N

for each j ∈ Iq , and we will set u = (u1 , u2 , . . . , uq ).

Similarly, notice that unj ∈ Vh can be represented as unj = (unj,k )k∈K for each j ∈ Iq . Also, we let u j = (unj )n∈I for each N j ∈ Iq , and u = (u 1 , u 2 , . . . , u q ). Definition 3.1. Define the inner product ·, · : Vh × Vh → R and the norm · 1 : Vh → R by

u, v = h∗

 k ∈K

uk vk ,

u 1 = h∗



|uk |,

∀u, v ∈ Vh .

(3.7)

k ∈K

The Euclidean norm induced by ·, · will be denoted by · 2 . Meanwhile, · ∞ will be the usual infinity norm in Vh .

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.6 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

6

Definition 3.2. For each sequence (vn )nN=0 ⊆ Vh and (k, n) ∈ K × I N −1 , we introduce the discrete linear operators:

⎧ vn+1 − vnk ⎪ ⎪ , ⎨ δt vnk = k ⎪ ⎪ ⎩ δ (2) vn = t

(1 )

δt vnk =

τ

vkn+1

k

− 2vnk + vkn−1 , 2

n t vk

=

μ

τ

vkn+1 − vkn−1 vkn+1

2τ

+ vnk

2

, (3.8)

.

It is well known that these difference operators provide consistent approximations to some differential operators under suitable analytical conditions. Moreover, we will set vˆ nk = μt vnk for briefness. Definition 3.3 (Ortigueira [61]). For any function f : R → R, any h > 0 and α of f at the point x is defined as ∞ 

h(α ) f (x) =

(α )

gk

f (x − kh),

α > −1, the fractional centered difference of order

∀x ∈ R,

(3.9)

k=−∞

whenever the double series at the right-hand side exists. Here, (α )

gk

=

(−1)k (α + 1)

, ( α2 − k + 1)( α2 + k + 1)

∀k ∈ Z.

(3.10)

Lemma 3.4 (Çelik and Duman [62]). If 1 < α ≤ 2 then the coefficients ( gk )k∞=−∞ satisfy: (α )

(α )

(a) g 0

(α )

(b) gk (c)

∞ 

> 0, (α ) = g −k ≤ 0 for all k = 0, and (α )

gk

(α )

= 0. As a consequence, it follows that g 0 = −

k=−∞

∞ 

(α )

gk .

k=−∞ k =0

Lemma 3.5 (Çelik and Duman [62]). Let f ∈ C 5 (R) and assume that all its derivatives up to order five are integrable. If 1 < α ≤ 2 then, for almost all x,

−

hα f (x) hα

=

dα f (x) d|x|α

+ O(h2 ).

(3.11)

Definition 3.6. For each sequence (vn )nN=0 ⊆ Vh , (k, n) ∈ K × I N −1 and i ∈ I p , define the linear operator Mi 1  (α ) n gk −l vk1 ,...,k

δx(αi ) vnk = − α h i

l =0

i

i −1 ,l,k i +1 ,...,k p

(3.12)

.

In light of Lemma 3.5, the operators (3.12) provide quadratically consistent approximations to the Riesz partial derivative of u with respect to xi of order α at the point (x1,k1 , . . . , xi −1,ki−1 , xi ,ki , xi +1,ki+1 , . . . , x p ,k p ) and time tn . Meanwhile, the fractional Laplacian will be estimated with a quadratic order of consistency, using the expression (α )

δx vnk =

p  (α ) δxi vnk .

(3.13)

i =1

The following results will be needed in our study. Lemma 3.7 (Macías-Díaz [40]). For each i ∈ I p and α ∈ (1, 2], there exists a unique positive self-adjoint (square-root) linear operator (α ) (α ) (α ) (α ) xi : Vh → Vh , such that −δxi u, v = xi u, xi v , for each u, v ∈ Vh . Lemma 3.8 (Macías-Díaz [63]). Let v ∈ Vh and i ∈ I p . Let α ∈ (1, 2] and define the constants

  p  −2 α (α ) (α ) g = 2g h∗  h , h

0

i

i =1

(α )

gh

(α )

= 2g 0 h∗

p  i =1

α h− . i

(3.14)

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.7 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

7

α (a) xi v 22 ≤ 2g 0 h∗ h−

v 22 . i (α )

(α )

(b) δxi

v 22

(α )

(α )

(α )

= xi xi v 22 .  2 (α ) (α ) α (c) δxi v 22 ≤ 4 g 0 h∗ h−

v 22 . i   (α )  2 (α ) (α ) (α )

δxi v 22 ≤ gh v 2 and

xi v 22 ≤ gh v 22 . (d) i ∈I p

i ∈I p

Definition 3.9. Suppose that G : Rq → R is a differentiable function, and let u j = (unj )n∈I ⊆ Vh , for each j ∈ Iq . For each N

j ∈ Iq and (k, n) ∈ K × I N −1 , we define the nonlinear operator

δt ,u j G nj,k (u) =

  ⎧ ˆ n1,k , . . . , uˆ nj−1,k , unj ,+k 1 , uˆ nj+1,k , . . . , uˆ nq,k − G j uˆ n1,k , . . . , uˆ nj−1,k , unj,k , uˆ nj+1,k , . . . , uˆ nq,k Gj u ⎪ ⎪ ⎪ ⎪ , ⎪ ⎪ ⎨ unj ,+k 1 − unj,k if unj ,+k 1 = unj,k , ⎪ ⎪ ⎪ ⎪ ∂ G (u (xk , tn+ 1 )) ⎪ ⎪ 2 ⎩ , ∂u j

if unj ,+k 1 = unj,k . (3.15)

Next, we extend dimensionally the operators of Definitions 3.2, 3.6 and 3.9 Definition 3.10. For each (k, n) ∈ K × I N , we let unk = (un1,k , un2,k , . . . , unq,k ). We define

 δt unk = δt un1,k , δt un2,k . . . , δt unq,k , ∀(k, n) ∈ K × I N −1 ,  (1 ) (1 ) (1 ) (1 ) δt unk = δt un1,k , δt un2,k . . . , δt unq,k , ∀(k, n) ∈ K × I N −1 ,  (2 ) (2 ) (2 ) (2 ) δt unk = δt un1,k , δt un2,k . . . , δt unq,k , ∀(k, n) ∈ K × I N −1 ,  μt unk = μt un1,k , μt un2,k , . . . , μt unq,k , ∀(k, n) ∈ K × I N −1 ,  (α ) (α ) (α ) (α ) δxi unk = δxi un1,k , δxi un2,k , . . . , δxi unq,k , ∀(k, n) ∈ K × I N , ∀i ∈ I p ,  δt ,u G nk (u) = δt ,u 1 G n1,k (u), δt ,u 2 G n2,k (u), . . . , δt ,uq G nq,k (u) , ∀(k, n) ∈ K × I N −1 .

(3.16) (3.17) (3.18) (3.19) (3.20) (3.21)

As expected, we will set

δx(α ) unk =

p  δx(αi ) unk ,

∀(k, n) ∈ K × I N .

(3.22)

i =1

In an analogous fashion, we let u nk = (un1,k , un2,k , . . . , unq,k ) for each (k, n) ∈ K × I N , and we can readily define the operators (3.16)–(3.22) for u as above. Using the hypothesis (H) and the notation introduced in the present section, the finite-difference scheme used in this work to solve (2.18) is the set of difference equations

μt δt(2) unk + γ δt unk = d⎧μt δx(α ) unk − δt ,u G nk (u),

∀(k, n) ∈ K × I N −2 , 0 u = φ( x ), ∀ k ∈ K , ⎪ k ⎪ ⎨ k1 uk = ψ(xk ), ∀k ∈ K, subject to ⎪ u2 = χ (xk ), ∀k ∈ K, ⎪ ⎩ nk uk = 0, ∀(k, n) ∈ ∂ K × I N .

(3.23)

Note that we chose to prescribe the initial data exactly by means of the functions φ, ψ, χ :  → Rq . Also, it is clear that this system is a fully explicit model, whence the existence and uniqueness of solutions is an immediate consequence. Indeed, for each (k, n) ∈ K × I N −1 , the difference equation of (3.23) can be rewritten as



unk +2 = unk +1 + unk − unk −1 + 2τ 2 dμt δx unk − δt ,u G nk (u) − γ δt unk (α )



(3.24)

The discrete model (3.23) is, thus, a four-step numerical scheme. For the sake of convenience, Fig. 1 provides the forwarddifference stencil of the discrete model in the case when p = q = 1.

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.8 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

8

Fig. 1. Forward-difference stencil for the approximation to the exact solution of the one dimensional form of (2.8) at the time tn , using the finite-difference scheme (3.23) when p = q = 1. The black circles represent the known approximations at the times tn−1 , tn and tn+1 , while the cross denotes the unknown approximation at the time tn+2 .

4. Numerical properties The properties of consistency, stability, boundedness and convergence of this scheme will be mathematically established in this section. In a first stage, we wish to prove the quadratic consistency of the scheme. To that end, we will consider the continuous operator L and the discrete operator L, which are respectively defined as

∂ 2 u (x, t ) ∂ u (x, t ) +γ − dα u (x, t ) + F (u (x, t )), ∀(x, t ) ∈ , ∂t ∂t2 (2 ) (α ) ∀(k, n) ∈ K × I N −2 . L u (unk ) = μt δt unk + γ δt unk − dμt δx unk + δt ,u G nk (u ),

Lu (u (x, t )) =

(4.1) (4.2)

Throughout, we will suppose that (H) is satisfied. Theorem 4.1 (Consistency). If u ∈ Cx,t () and F ∈ L ∞ (Rq ) then there exists a constant C which is independent of h and τ , such that 5, 4

Lu (u (xk , tn+ 1 )) − 2

L u (unk ) ∞

≤ C (τ 2 + h 22 ), for each (k, n) ∈ K × I N −2 .

Proof. Note that the smoothness of u and the essential boundedness of F imply that u j ∈ Cx,t () and F j ∈ L ∞ (Rq ), for each j ∈ Iq . The consistency of the discrete operators in Definitions 3.2, 3.6 and 3.9 along with Taylor’s theorem guarantee 5, 4

j

j

j

j

τ , such that   2   ∂ u j (xk , tn+ 1 )   j (2 ) 2 − μ δ u ( x , t ) ∀(k, n) ∈ K × I N −2 , ∀ j ∈ Iq ,  ≤ C1τ 2,  t n j k t   ∂t2     ∂ u j (xk , tn+ 1 )   j 2 − δt u j (xk , tn ) ≤ C 2 τ 2 , ∀(k, n) ∈ K × I N −1 , ∀ j ∈ Iq ,    ∂t   α   ∂ u j (xk , tn+ 1 )   j (α ) 2 − μ δ u ( x , t ) ∀(k, n) ∈ K × I N −1 , ∀(i , j ) ∈ I p × Iq ,  ≤ C 3,i (h2i + τ 2 ),  t n j k x i   ∂|xi |α     j ∀(k, n) ∈ K × I N −1 , ∀ j ∈ Iq .  F j (xk , tn+ 1 ) − δt ,u j G nk (u ) ≤ C 4 τ 2 ,

now that there exist constants C 1 , C 2 , C 3,i and C 4 which are independent of h and

(4.3)

(4.4)

(4.5) (4.6)

2

j

j

Let C l = max{C l : j ∈ Iq } for each l = 1, 2, 4, and define C 3 = max{C 3,i : (i , j ) ∈ Iq × Iq }. Clearly, these constants are independent of h and τ . Moreover, using the triangle inequality, it follows that

    Lu (u (xk , tn+ 1 )) − L u (unk ) 2

 2     ∂ u (xk , tn+ 1 )   ∂ u (xk , tn+ 1 )    (2 ) n  n 2 2 ≤ − μ δ u + γ − δ u    t t t k k     ∞ ∂t ∂t2 ∞ ∞        (α )  + d α u (xk , tn+ 1 ) − μt δx unk  +  F (u (xk , tn+ 1 )) − δt ,u G nk (u ) 2

∞

≤ C 1 τ 2 + γ C 2 τ 2 + dC 3 ( h 22 + τ 2 ) + C 4 τ 2 ,

2

(4.7) ∞

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.9 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

9

for each (k, n) ∈ K × I N −2 . The conclusion of the theorem readily follows with C = max{C 1 , γ C 2 , dC 3 , C 4 }.

2

Lemma 4.2 (Macías-Díaz [40]). Let G : Rq → Rq be such that G ∈ C 2 (Rq ) and D 2 G ∈ L ∞ (R), and let (unj )nN=0 , (vnj )nN=0 and (Rnj )nN=0

be sequences in Vh , for each j ∈ Iq . Let enj = vnj − unj and  G nj = δ v j ,t G n (v) − δu j ,t G n (u), for each n ∈ I N −1 and j ∈ Iq . There exist constants C 1 , C 2 , C 3 ∈ R+ which depend only on G, such that, for each j ∈ Iq

 G nj 22 ≤ C 1 ( enj +1 22 + enj 22 ), 2τ

m  m m       n n

Rnj 22 + C 2 e0j 22 + C 3 τ

δt enj 22 , R j − G j , δt enj  ≤ 2τ n =1

n =0

mτ 2

m 

∀n ∈ I N −1 ,

(4.8)

∀m ∈ I N −1 ,

(4.9)

∀m ∈ I N −1 .

(4.10)

n =0

 G nj 22 ≤ 4C 1 T 2 e0j 22 + 4C 1 T 3 τ

n =1

m 

δt enj 22 ,

n =0

Lemma 4.3 (Macías-Díaz [40]). Let G, (unj )nN=0 , (vnj )nN=0 and (Rnj )nN=0 be as in Lemma 4.2. Let enj = vnj − unj and  G nj = δ v j ,t G n (v) −

δu j ,t G n (u), for each n ∈ I N −1 and j ∈ Iq . If

μt δt(2) enj + γ δt enj − dμt δx(α ) enj +  G nj = Rnj ,

∀n ∈ I N −1 , ∀ j ∈ Iq ,

(4.11)

then for each m ∈ I N −1 and j ∈ Iq , p    (α ) (α )

xi e1j 22 + xi emj +1 22

τ 2 δt(2) emj +1 22 ≤ 160C 1 T 2 e0j 22 + 20μt δt e0j 22 + 5p τ 2 gh(α ) d2

i =1

+40T τ

m 

Rnj 22

2

2

+ 20(8C 1 T + γ ) T τ

n =1

m 

(4.12)

δt enj 22 .

n =0

The following discrete version of Gronwall’s inequality will be needed in the sequel. Lemma 4.4 (Pen-Yu [64]). Let (ωn )nN=0 and (ρ n )nN=0 be finite sequences of nonnegative mesh functions, and suppose that there exists C ≥ 0 such that

ωk ≤ ρ k + C τ

k −1 

ωk , ∀k ∈ IN −1 .

(4.13)

n =0

Then ωn ≤ ρ n e Cnτ for each n ∈ I N .

2

Theorem 4.5 (Stability). Let G : Rq → Rq be such that G ∈ C 2 (Rq ) and D 2 G ∈ L ∞ (R), and suppose that u and v are solutions of (α ) (3.23) corresponding to the initial conditions (φu , ψu , χu ) and (φ v , ψ v , χ v ), respectively. Let e = u − v, and assume that 52 p τ 2 gh < + 1. Then there exists a constant C 0 ∈ R which is independent of τ , h, u and v, and a constant 0 < η0 < 1 independent of u and v, such that

1 2

δt enj 22

Proof. Let

+ (1 − η0 )d

p 

(α ) n 2

xi e j 2 ≤ C 0

e0j 22

0 2 t δt e j 2

+μ

i =1

+

p 



(α ) 1 2

xi e j 2 ,

∀n ∈ I N −1 , ∀ j ∈ Iq .

(4.14)

i =1

η0 satisfy 52 p τ 2 gh(α ) < η0 < 1. Note that the following discrete problem is satisfied:

(α ) μt δt(2) enk + γ δt enk − dμt ⎧ δx enk + δt ,u G nk (u) − δt , v G nk (v) = 0, 0

e = φu (xk ) − φ v (xk ), ⎪ ⎪ ⎨ k1 ek = ψu (xk ) − ψ v (xk ), subject to ⎪ ek2 = χu (xk ) − χ v (xk ), ⎪ ⎩ n ek = 0,

∀(k, n) ∈ K × I N −2 , ∀k ∈ K, ∀k ∈ K, ∀k ∈ K, ∀(k, n) ∈ ∂ K × I N .

(4.15)

Following the nomenclature introduced at the beginning of this section, we identify the left-hand side of the difference equation of (4.15) as L u (unk ) − L v (vnk ). Moreover, for the sake of convenience, we let  G nk = δt ,u G nk (u) − δt , v G nk (v), for each

(k, n) ∈ K × I N −1 . For each j ∈ Iq and n ∈ I N −2 , the following identities are satisfied:

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.10 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

10

1

τ2

2 1

4

μt δt(2) enj , δt enj = δt μt δt enj −1 22 − (α )

(α )

−μt δxi enj , δt enj = δt xi enj 22 ,

δt δt(2) enj 22 ,

(4.16)

∀i ∈ I p .

2

(4.17)

Next, we calculate the inner product of δt enj with both sides of the difference equation L u (unj ) − L v (vnk ) = 0, substitute then the identities above, and we obtain next the sum of the resulting identity for all n ∈ Im . Multiply then by 2τ on both sides, apply Lemma 4.2 with Rn = 0 and simplify algebraically to obtain that, for each m ∈ I N −1 ,

1 2

δt emj +1 22 + d

p 

xi emj +1 22 ≤ μt δt emj 22 + d (α )

i =1 p 

(α )

xi e1j 22 +

i =1

5 2

xi emj +1 22 (α )

i =1

≤ μt δt e0j 22 + d ≤ ρj +

p 

(α )

p τ 2 gh d2

p 

τ

2

2

δt emj +1 22 + 2τ (2 )

m      n G j , δt enj 

(4.18)

n =1

(xαi ) emj +1 22 + C 5 τ

m 

δt enj 22 ≤ ρ j + dη0

n =0

i =1

p 

(xαi ) emj +1 22 + C 5 τ

m 

ωnj .

n =0

i =1

Here, C 4 = max{C 2 + 80C 1 T , 11}, C 5 = 2C 3 + 20(8C 1 T + γ ) T and 2

2

ρ j = C 4 e0j 22 + μt δt e0j 22 + d

p 

2



(α ) 1 2

xi e j 2 ,

(4.19)

i =1

1

ωnj = δt enj 22 + (1 − η0 )d 2

p 

(α )

xi enj 22 ,

∀n ∈ I N −1 .

(4.20)

i =1

Subtracting the second term on the right-hand side of (4.18), we note that the hypotheses of Lemma 4.4 are readily satisfied with C = C 5 and ρ k = ρ for each k ∈ I N −2 . The conclusion of this theorem follows with C 0 = C 4 e C 5 T . 2 The following result is readily obtained using a proof similar to that of Theorem 4.5. Corollary 4.6 (Boundedness). Let G : Rq → Rq satisfy G ∈ C 2 (Rq ) and D 2 G ∈ L ∞ (R). If u is the solution of (3.23) then there exists a constant C ∈ R+ such that unk ∞ ≤ C , for each (k, n) ∈ K × I N . 2 Theorem 4.7 (Convergence). Let u ∈ Cx,t () be a solution of (2.18), and suppose that G ∈ C 2 (R) and D 2 G ∈ L ∞ (R). If 5, 4

1 and (3.23) has exact initial data then the numerical solution converges to that of the continuous problem with order in the Euclidean norm. Proof. The proof is similar to that of Theorem 4.5. Fix for each (k, n) ∈ K × I N , and define

 = n j

unj

− unj .

(α ) 5 p 2 gh d < 2 2 O( + h 22 )

τ τ

η0 as in that theorem, let Rnk be the local truncation error at (xk , tn ),

Then the following system is satisfied: n

) n μt δt(2)  nk + γ δt  nk − dμt δx(α  k + δt ,u G nk (u) − δt ,u G nk (u ) = Rk , ∀(k, n) ∈ K × I N −2 ,  k0 =  k1 =  k2 = 0, ∀k ∈ K, subject to  nk = 0, ∀(k, n) ∈ ∂ K × I N .

(4.21)

Mimicking the proof of Theorem 4.5, we let  G nk = δt ,u G nk (u) − δt ,u G nk (u ), for each (k, n) ∈ K × I N −2 . Proceeding as in that proof of the previous result, we reach the inequality

1 2

+1 2

δt m

2 j



+d

p 

( α ) m +1 2

xi  j

2 ≤ ρ

m +1

+

i =1

5 2

pτ

2 (α ) gh d

p 

+1 2

xi  m

2 + C 5 τ j (α )

m 

ωn , ∀k ∈ IN −1 ,

(4.22)

n =0

i =1

where C 4 = max{C 2 + 80C 0 T 2 , 11, 20T + 2}, the constant C 5 is as in the proof of Theorem 4.5 and

m j

ρ = C 4 

0 2 j 2

0 2 t δt j 2

+μ



+d

p 

(α ) 1 2

xi  j 2 + τ

i =1 m j

ω =

1 2

2

δt m j 2



+ (1 − η0 )d

p  i =1

2

(xαi )  m j 2 ,

m −1 



Rnj 22

,

∀m ∈ I N −1 ,

(4.23)

∀m ∈ I N −1 .

(4.24)

n =0

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.11 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

11

Move the second term on the right-hand side of (4.22) to the left-hand side and apply Lemma 4.4. Use then the initial data of (4.21), and use Theorem 4.1 to see that there exists a constant A 1 independent of τ and h, such that

1 2

δ 

n 2 j 2

≤

1 2

δt nj 22



+ (1 − η0 )d

p 

(α )

xi  nj 22 ≤ A 1 τ

m −1 

Rnj 22 ≤ A 2 (τ 2 + h 22 )2 ,

∀n ∈ I N −1 , ∀ j ∈ Iq .

n =0

i =1

(4.25) Here, A 2 = A 1 C T , and C is the constant provided by Theorem 4.1. Multiply both ends of (4.25) by 2 and use the reversed √ triangle inequality to obtain that  nj +1 2 −  nj 2 ≤ 2 A 2 τ (τ 2 + h 22 ). Next, take the sum on both sides of this inequality for n between 0 and m − 1, the formula for telescoping sums and the data at n = 0 to obtain

 m j 2 ≤ where A 0 =

√



2 A2τ

m −1 

(τ 2 + h 22 ) ≤ A 0 (τ 2 + h 22 ),

∀m ∈ I N , ∀ j ∈ Iq ,

(4.26)

n =0

2 A 2 T . The conclusion readily follows now.

2

Remarks 4.8. Before closing this section, it is important to point out that the stability and the convergence of the finitedifference scheme has been established using a discrete and fractional form of the well-known energy method. This approach has been followed by the authors in some previous papers in order to analyze some Hamiltonian systems consisting of a single hyperbolic partial differential equation with fractional diffusion [33,40,55]. However, most of the hyperbolic models in which complex patterns appear, are non-Hamiltonian regimes consisting of two or more partial differential equations. In that sense, the previous efforts by this author are useless in the investigation of such systems. We have proved rigorously that the present methodology is numerically efficient, and it has the advantage of being applicable to a wide class of systems consisting of various hyperbolic fractional partial differential equations. 5. Computer implementation In this section, we will describe an efficient computational implementation of the finite-difference method (3.23). To that end, we will focus our attention on the explicit equations (3.24), and we will set p = 3 and q = 2. In order to simplify the notation, we will convey that x = x1 , y = x2 and z = x3 , and that u = u 1 and v = u 2 . Moreover, suppose that τ1 = τ2 = 1, γ1 = γ2 = γ ∈ R+ ∪ {0}, d1 = d2 = d ∈ R+ ∪ {0} and α1 = α2 = α ∈ (1, 2]. Using these conventions and the hypothesis (H), the mathematical model (2.8) can be rewritten as

∂ 2 u (x, t ) ∂ u (x, t ) +γ = dα u (x, t ) − F 1 (u (x, t ), v (x, t )), ∀(x, t ) ∈ , 2 ∂t ∂t 2 ∂ v (x, t ) ∂ v (x, t ) +γ = dα v (x, t ) − F 2 (u (x, t ), v (x, t )), ∀(x, t ) ∈ , ∂t ∂t2 ⎧ u (x, 0) = φ1 (x) v (x, 0) = φ2 (x), ∀x ∈ B , ⎪ ⎨ ∂ u (x, 0) ∂ v (x, 0) such that = ψ1 (x), = ψ2 (x), ∀x ∈ B , ⎪ ∂t ∂t ⎩ u (x, t ) = v (x, t ) = 0, ∀(x, t ) ∈ ∂ B × [0, T ],

(5.1)

where

∂ G 1 (u (x, t ), v (x, t )) , ∂u ∂ G 2 (u (x, t ), v (x, t )) , F 2 (u (x, t ), v (x, t )) = ∂v

F 1 (u (x, t ), v (x, t )) =

∀(x, t ) ∈ ,

(5.2)

∀(x, t ) ∈ .

(5.3)

Let also k = k1 , l = k2 and m = k3 , let j = (k, l, m) and agree that xk = x1,k , yl = x2,l and zm = x3,m . Under the new conventions, if ( j , n) ∈ K × I N −2 then the recursive formulas (3.24) can be rewritten equivalently as

⎧   G 1 (unj +1 , μt vnj ) − G 1 (unj , μt vnj ) ⎪ ⎪ ( α ) n +1 (α ) n n +2 n +1 n −1 n ⎪ u = u + u − u + r u + δ u r − r3 unj +1 − unj , δ − ⎪ 1 2 x x j j j j j j ⎪ n + 1 n ⎨ u −u j

j

⎪   G 2 (μt unj , vnj +1 ) − G 2 (μt unj , vnj ) ⎪ ⎪ ( α ) n +1 (α ) n n +2 n +1 n −1 n +1 n n ⎪ − r v − v , ⎪ 3 j j ⎩ v j = v j + v j − v j + r 1 δx v j + δx v j − r 2 vnj +1 − vnj (5.4) where r1 = τ 2 d, r2 = 2τ 2 and r3 = 2τ γ . Meanwhile, the initial-boundary data will assume the form

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.12 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

12

⎧ 0 u = φ u (xk ), vk0 = φ v (xk ), ⎪ ⎪ ⎨ k1 uk = ψ u (xk ), vk1 = ψ v (xk ), ⎪ u2 = χ u (xk ), vk2 = χ v (xk ), ⎪ ⎩ nk uk = vnk = 0,

∀k ∈ K, ∀k ∈ K, ∀k ∈ K, ∀(k, n) ∈ ∂ K × I N .

(5.5)

Before we propose the computational implementation of the system (5.4), it is important to point out that its extension to account for other values of q is straightforward. In what follows, for any α ∈ (1, 2] and n ∈ N , we consider the real matrices of sizes (n + 1) × (n + 1) defined by

⎛

(α )

(α )

g0

⎜ g (α ) ⎜ 1 ⎜ (α ) g H nα = ⎜ ⎜ 2 ⎜ .. ⎝ .

(α )

g −2

g −3

g0

g −1

g −2

g1

g0

g −1

(α )

(α )

(α )

(α )

(α )

.. .

(α )

.. .

(α )

g n −1

g n −2

(α )

··· ··· ··· .. . ···

(α )

(α )

.. .

gn

(α )

g −1

(α )

g n −3

g −n

⎞

(α ) g 1−n ⎟ ⎟

⎟ ⎟ ⎟ ⎠

(α ) g 2−n ⎟ .

.. .

(5.6)

(α )

g0

We will enumerate the rows and columns of matrices beginning from 0, and not from 1. It is easy to see then that the (α ) component of H nα at the ith row and jth column is given by [ H nα ]i , j = g i − j , for all (i , j ) ∈ I n × I n . Moreover, the properties of the coefficients ( gk )k∞=−∞ summarized in Lemma 3.4 assure that the matrix H nα is Hermitian. Let w ∈ R M 1 +1 × R M 2 +1 × R M 3 +1 . If k ∈ I M 1 then we agree that wk,·,· ∈ R M 2 +1 × R M 3 +1 represent the matrix (α )

⎛

wk,0,0 wk,1,0 wk,2,0

wk,0,1 wk,1,1 wk,2,1

wk,0,2 wk,1,2 wk,2,2

.. .

.. .

.. .

wk, M 2 ,0

wk, M 2 ,1

wk, M 2 ,2

⎜ ⎜ ⎜ wk,·,· = ⎜ ⎜ ⎝

··· ··· ··· .. . ···

wk,0, M 3 wk,1, M 3 wk,2, M 3

.. .

⎞

⎟ ⎟ ⎟ ⎟. ⎟ ⎠

(5.7)

wk, M 2 , M 3

In other words, wk,·,· is obtained from w by fixing the first component equal to k, and arranging naturally the entries as an ( M 2 + 1) × ( M 3 + 1) matrix. In an entirely similar fashion, if l ∈ I M 2 and m ∈ I M 3 then w·,l,· ∈ R M 1 +1 × R M 3 +1 and w·,·,m ∈ R M 1 +1 × R M 2 +1 are the matrices of sizes ( M 1 + 1) × ( M 3 + 1) and ( M 1 + 1) × ( M 2 + 1) given respectively by

⎛

w0,l,0 w1,l,0 w2,l,0

w0,l,1 w1,l,1 w2,l,1

w0,l,2 w1,l,2 w2,l,2

.. .

.. .

.. .

w M 1 ,l,0

w M 1 ,l,1

w M 1 ,l,2

⎜ ⎜ ⎜ w·,l,· = ⎜ ⎜ ⎝ and

⎛ ⎜ ⎜ ⎜ w·,·,m = ⎜ ⎜ ⎝

··· ··· ··· .. . ···

w0,0,m w1,0,m w2,0,m

w0,1,m w1,1,m w2,1,m

w0,2,m w1,2,m w2,2,m

.. .

.. .

.. .

w M 1 ,0,m

w M 1 ,1,m

w M 1 ,2,m

⎞

w0,l, M 3 w1,l, M 3 w2,l, M 3

⎟ ⎟ ⎟ ⎟ ⎟ ⎠

.. .

(5.8)

w M 1 ,l, M 3

··· ··· ··· .. . ···

w0, M 2 ,m w1, M 2 ,m w2, M 2 ,m

.. .

⎞ ⎟ ⎟ ⎟ ⎟. ⎟ ⎠

(5.9)

w M 1 , M 2 ,m

We will describe next the implementation of the calculations involving the terms with fractional centered differences. To that end, we let w be any of u or v as before, and recall that (3.13) holds. The first term on the right-hand side of that expression can be equivalently rewritten as M1 M1 1  (α ) n 1  α α α gk− j w j ,l,m = − α [ H M 1 ]kj [wn·,·,m ] jl = −h− 1 [ H M 1 w·,·,m ]k,l . h1 1

δx(α ) wnk,l,m = − α h

j =0

(5.10)

j =0

Obviously, the product inside of the parenthesis at the right-hand side of this expression is the usual product of matrices. (α ) (α ) n α α −α α n n Similarly, it is easy to check that δ y wnk,l,m = −h− 2 [ H M 2 wk,·,· ]l,m and δ z wk,l,m = −h 3 [wk,·,· H M 3 ]l,m . Combining these expressions and representing the discrete fractional Laplacian operator simply by δ (α ) , it follows that α α −α −α n α n n α δ (α ) wnk,l,m = −h− 1 [ H M 1 w·,·,m ]k,l − h 2 [ H M 2 wk,·,· ]l,m − h 3 [wk,·,· H M 3 ]l,m ,

for each (k, l, m) ∈ K.

(5.11)

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.13 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

13

Algorithm 1: Iterative algorithm to solve the discrete problem (3.23). 1 Initialize parameters; 2 Define matrix H ; 3 Set uk0,l,m ← φku,l,m ; vk0,l,m ← φkv,l,m ; 4 Set uk1,l,m ← ψku,l,m ; vk1,l,m ← ψkv,l,m ; 5 Set uk2,l,m ← χku,l,m ; vk2,l,m ← χkv,l,m ; 6 for j = 0 : M do 7 for w = u , v do w 8 Set R ·,·, = Hw1·,·, j ; j 1 Set S w j ,·,· = Hw j ,·,· ;

9

Set T w = w1j ,·,· H ; j ,·,·

10 11

end

12 end 13 for n = 1 : N − 2 do

G 1 (unk,+l,1m , μt vnk,l,m ) − G 1 (unk,l,m , μt vnk,l,m )

14

Set W k,l,m ←

15 16

for j = 0 : M do for w = u , v do +1 w Set r·,·, = Hwn·,·, ; j j

17

unk,+l,1m

− unk,l,m

; Z k,l,m ←

G 2 (μt unk,l,m , vnk,+l,1m ) − G 2 (μt unk,l,m , vnk,l,m ) vnk,+l,1m − vnk,l,m

+1 Set s w = Hwnj ,·,· ; j ,·,·

18

n+1 Set t w j ,·,· = w j ,·,· H ;

19

end

20 21

end

22

Calculate un+2 = un+1 + un − un−1 − r3 (un+1 − un ) − r2 W − r u − su − t u − R u − S u − T u ;

23 24 25

;

Calculate vn+2 = vn+1 + vn − vn−1 − r3 (vn+1 − vn ) − r2 Z − r v − s v − t v − R v − S v − T v ; Set R u ← r u ; S u ← su ; T u ← t u ; Set R v ← r v ; S v ← s v ; T v ← t v ;

26 end

α α −α α −α α We will require the matrices H x = r1 h− 1 H M 1 , H y = r 1 h 2 H M 2 and H z = r 1 h 3 H M 3 . Using this nomenclature, the finitedifference system (5.4) way be rewritten equivalently as

⎧  G 1 (ukn,+l,1m , μt vnk,l,m ) − G 1 (unk,l,m , μt vnk,l,m ) ⎪ ⎪ n +2 n +1 n −1 n ⎪ u = u + u − u − r − r3 ukn,+l,1m − unk,l,m , ⎪ 2 k,l,m k,l,m k,l,m k,l,m n + 1 ⎪ n ⎪ uk,l,m − uk,l,m ⎪ ⎪ ⎪   ⎪ ⎪ n +1 n +1 +1 n n n ⎪ ⎨ − [ H x un·,·, ] + [ H u ] + [ u H ] H u ] + [ H u ] + [ u H ] − [ , y z x y z k , l l , m l , m k , l l , m l , m m ·,·,m k,·,· k,·,· k,·,· k,·,· ⎪  G 1 (μt unk,l,m , vkn,+l,1m ) − G 1 (μt unk,l,m , vnk,l,m ) ⎪ ⎪ n +2 n +1 n −1 n +1 n n ⎪ v = v + v − v − r − r v − v ⎪ 2 3 k,l,m k,l,m , k,l,m k,l,m k,l,m ⎪ k,l,m ⎪ vkn,+l,1m − vnk,l,m ⎪ ⎪ ⎪   ⎪ ⎪ n +1 n +1 +1 n n n ⎩ − [ H x vn·,·, m ]k,l + [ H y vk,·,· ]l,m + [vk,·,· H z ]l,m − [ H x v·,·,m ]k,l + [ H y vk,·,· ]l,m + [vk,·,· H z ]l,m ,

(5.12)

for each (k, l, m, n) ∈ K × I N −2 . In what follows, we will describe the computer implementation when M = M 1 = M 2 = M 3 , a = a1 = a2 = a3 and b = b1 = b2 = b3 . In such case, we observe that h = h1 = h2 = h3 and H = H x = H y = H z . Let n ∈ I N −2 and w = u , v, and define the three-dimensional real arrays R w , S w and T w of sizes ( M + 1) × ( M + 1) × ( M + 1), by

R kw,l,m = [ Hwn·,·,m ]k,l ,

S kw,l,m = [ Hwnk,·,· ]l,m ,

T kw,l,m = [wnk,·,· H ]l,m ,

(5.13)

for all (k, l, m) ∈ K. For simplification purposes, we are obviating here the dependence of R w , S w and T w on k. From these set of identities, it readily follows that for each (k, l, m) ∈ K, w n R ·,·, m = Hw·,·,m ,

S kw,·,· = Hwnk,·,· ,

T kw,·,· = wnk,·,· H .

(5.14)

Under these circumstances, Algorithm 1 provides the description of our computational implementation of (5.12). Remarks 5.1. Some important remarks on the computer implementation of Algorithm 1 must be mentioned. 1. Firstly, it is important to point out that our implementation was carried out having Fortran 95 in mind. However, the approach is also valid for other computer languages, like C++ or Matlab. We opted to use Fortran in view that it is a free computer language, and our experience tells us that it is usually faster than other languages.

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.14 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

14

2. Various code lines can be obviously executed in parallel. One example is the set of instructions between lines 15–20 of Algorithm 1. Also, lines 22 and 23 can be realized in parallel, as well as all the instruction in lines 24 and 25. All these facts have been exploited in our computer implementation of Algorithm 1. 3. The algorithm was coded in Gfortran on a Linux Mint 18.3 “Sylvia” MATE distribution, and the use of Fortran commands for matrix multiplication has been crucial to speed up our implementation. In turn, the parallelization was carried out naturally using the package OpenMP. (α ) 4. In the calculation of the coefficients ( gk )k∞=−∞ , we employed the following recursive formula:

⎧ (α + 1) (α ) ⎪ ⎪ ⎨ g 0 = (α /2 + 1)2 ,   ⎪ (α ) α+1 (α ) ⎪ ⎩ g k +1 = 1 − g , α /2 + k + 1 k

(5.15)

∀k ∈ N ∪ {0}.

6. Illustrative simulations The present section is devoted to provide some illustrative computer simulations of the methods proposed in the previous sections. We will consider firstly the case of the non-variational scheme (3.23), and provide some three-dimensional simulations to the investigation of Turing patterns in chemical systems of inhibitor-activator substances. More concretely, consider the system

∂ 2 u (x, t ) ∂ u (x, t ) + = u (x, t ) − av (x, t ) + bu (x, t ) v (x, t ) − [u (x, t )]3 + du ∇ α1 u (x, t ), ∀(x, t ) ∈ , ∂t ∂t2 ∂ 2 v (x, t ) ∂ v (x, t ) τv + = u (x, t ) − cv (x, t ) + d v ∇ α2 v (x, t ), ∀(x, t ) ∈ , (6.1) ∂t2 ! ∂t u (x, 0) = φ u (x), v (x, 0) = φ v (x), ∀x ∈ B , ∂v ∂u such that (x, 0) = ψ u (x), (x, 0) = ψ v (x), ∀x ∈ B . ∂t ∂t It is clear that this is a particular form of (2.8), using q = 2, u = u 1 , v = u 2 , τu = τ1 , τ v = τ2 , γ1 = γ2 = 1, φ u = φ1 , φ v = φ2 ,

τu

ψ u = ψ1 and ψ v = ψ2 . Meanwhile, the functions F 1 and F 2 are provided by (2.12) and (2.13), respectively. Clearly, the functions G 1 , G 2 : R2 → R defined below satisfy the hypothesis (H): G 1 (u , v ) = 12 u 2 − auv + 12 bu 2 v − 14 u 4 , G 2 (u , v ) = uv −

1 cv 2 , 2

∀u , v ∈ R,

∀u , v ∈ R.

(6.2) (6.3)

In the following examples, we concentrate our attention on the solutions of the system (6.1). Our examples are motivated by a set of results reported in [44] for the two-dimensional version of the mathematical model considered here. Briefly, the authors of that work found out that, for any value of α1 , α2 ∈ (1, 2], Turing patterns appear in the system when the following parameter conditions are satisfied:

1 < c < a < aT , 1 (d 4

d > 1,

b<



c (a − c ).

(6.4)

2 −1

Here, a T = + c ) d . Motivated by this fact, we will fix the parameter values recorded in Table 1. We will only let b take on various values in R+ . Example 6.1. Consider the system (6.1) with the parameter values of Table 1. In addition, we will use the parameter value b = 0.5. As initial profiles, we chose samples of a uniformly distributed random variable on the interval [−0.03, 0.03], and used zero initial velocities. More precisely, the approximations at the time t = 0 will be samples of a random variable with uniform distribution on [−0.03, 0.03]. This task can be easily achieved by using a suitable form of the ran command in Fortran 95, in order to define the initial approximations uk0 , for each k ∈ K. To account for the zero initial velocities, we will

let uk2 = uk1 = uk0 , for each k ∈ K. Under those circumstances, Fig. 2 shows snapshots of the approximate solutions of the variable u of (2.8) versus x and y for α = 1.5. The results were obtained using the scheme (3.23) with the computational parameters in Table 1. The times (a) t = 0, (b) t = 100, (c) t = 400 and (d) t = 2500 were used in this example, and the graphs have been normalized with respect to the absolute maximum of the solution at each time. It is clear that the results are in good qualitative agreement with those obtained in [44]. Indeed, notice that the snapshots depict the evolution from an almost spatially homogeneous state of the system to the stationary spatially non-homogeneous state. Turing structures are obtained as a result. More precisely, in this case we witness the coexistence of stripes with spheres, as predicted in [44]. It is worth pointing out that we have conducted more experiments. The results are not shown in this work in order to avoid redundancy, but they are in agreement with the literature. 2 Our next example considers the three dimensional version of (6.1).

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.15 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

15

Table 1 Set of fixed model and computational parameters employed in the simulations of Examples 6.1 and 6.2. Parameter

p

Example 6.1 Example 6.2

2 3

a 7.45

c 5

du 1

dv 20

τu 1

τv 1

α1 1.5

α2



T

h

τ

1.5

[0, 200]3 [0, 100]3

2500 2000

1

0.01

Fig. 2. (Color online.) Snapshots of the approximate solutions of the variable u of (6.1) versus (x, y ), using the parameters of Table 1 and b = 1.5. The approximations at the time t = 0 were samples of a random variable with uniform distribution on [−0.03, 0.03], and we considered zero initial velocities. The times (a) t = 0, (b) t = 100, (c) t = 400 and (d) t = 2500 were used in these simulations. The graphs were normalized with respect to the absolute maximum of the solution at each particular time.

Example 6.2. Consider the problem (6.1) in three spatial dimensions, together with the set of model and computational parameters in Table 1. Additionally, we let b = 1.5. Fig. 3 shows snapshots of the approximate solution obtained using the scheme (3.23) at various times, namely, (a) t = 50, (b) t = 100, (c) t = 250 and (d) t = 2000. The results depict the presence of complex patterns in this system. The results were obtained using a parallel implementation of our code in Fortran 95, and the visualizations were obtained using standard routines in Matlab. In particular, we employed the standard command isosurface with different values of the isovalue parameter to that end. More precisely, we used values of the parameter isovalue equal to 0.01. Fig. 4 shows some x-, y- and z-cross sections of the solutions at each of those times. The solutions exhibit some two-dimensional complex patterns on each of the sides of the cubes. As in the previous examples, the approximations at the time t = 0 were samples of a random variable with uniform distribution on [−0.03, 0.03], and we imposed zero initial velocities. 2 Remarks 6.3. As we mentioned at the end of Section 3, the models in Example 2.3 satisfy the hypothesis (H). This fact was observed at the beginning of the present section for the model of Example 2.3(c), and we will show next that the other two examples satisfy this analytical condition.

• In the case of Example 2.3(a), a function satisfying the hypothesis is G 1 (u ) = 1 − cos(u ), for each u ∈ R.

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.16 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

16

Fig. 3. (Color online.) Approximation to the solutions of (6.1), using the finite-difference scheme (3.23) with the parameters of Table 1 and b = 1.5. The approximations at the time t = 0 were samples of a random variable with uniform distribution on [−0.03, 0.03], and we considered zero initial velocities. Various times were considered, namely, (a) t = 50, (b) t = 100, (c) t = 250 and (d) t = 2000. The plotting subroutines were based on the function isosurface with values of the isovalue parameter equal to 0.01.

• For Example 2.3(b) with the scaled reaction functions (2.10)-(2.11) the functions G 1 (u 1 , u 2 ) =

1 S

u 22 ln(u 2 + Su 1 ) − u 1 u 2 +

R 2

G 2 (u 1 , u 2 ) = Su 1 u 2 − S 2 u 21 ln(u 2 + Su 1 ) −

u 21 − Q 2

R 2S

u 22 ,

u 31 ,

∀u 1 , u 2 ∈ R,

∀u 1 , u 2 ∈ R,

(6.5) (6.6)

satisfy the hypothesis (H). Also, we must mention that the simulations of Example 6.2 show the appearance of complex patterns. However, it is important to point out that the analytical prediction of patterns in the three-dimensional model (6.1) is still an open task of research. In that sense, the present computational model can be a useful tool in the resolution of such problem. 7. Conclusion In this work, we investigated numerically a system of hyperbolic partial differential equations with fractional diffusion and coupled reaction terms. The mathematical model has various applications, depending on the form of the reaction functions. Motivated by these facts, we proposed an explicit finite-difference models to approximate the solutions of the continuous model. The discrete model is a non-variational scheme for which the properties of consistency, stability, boundedness and quadratic convergence are rigorously proved. To that end, an analytical constraint is imposed on the reaction terms of the mathematical model, and a discrete fractional form of the energy method is employed in order to establish rigorously the stability and the quadratic convergence of the scheme. The discretization is based on the use of fractional centered differences, and the computational implementation is carried out using parallel computing in Fortran 95. To that end, a convenient vector reformulation of the numerical method is proposed. Some illustrative applications of our methodology were presented in this work. Indeed, the finite-difference scheme was employed to solve some nonlinear systems which present Turing patterns in the two-dimensional scenario. Our simulations reflect this fact. Moreover, we used a

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.17 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

17

Fig. 4. (Color online.) Approximation to some x-, y- and z-cross sections of the solutions of (6.1), using the finite-difference scheme (3.23) with the parameters of Table 1 and b = 1.5. The approximations at the time t = 0 were samples of a random variable with uniform distribution on [−0.03, 0.03], and we considered zero initial velocities. Various times were considered, namely, (a) t = 50, (b) t = 100, (c) t = 250 and (d) t = 2000.

three-dimensional implementation of our scheme to obtain Turing patterns in some three-dimensional systems. It is well known that the resolution of fractional systems in three-dimensional scenarios is computationally highly demanding. However, our computer implementation was able to carry out this task, and exhibited the presence of complex patterns in the media considered herein. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement The author wishes to thank the anonymous reviewers for their comments and criticisms. All their suggestions were literally followed, and they resulted in a substantial improvement in the overall quality of this work. The present work reports on a set of final results of the research project “Conservative methods for fractional hyperbolic systems: analysis and applications”, funded by the National Council of Science and Technology, Mexico through grant A1-S-45928. References [1] V. Dufiet, J. Boissonade, Dynamics of Turing pattern monolayers close to onset, Phys. Rev. E 53 (5) (1996) 4883. [2] A. De Wit, Spatial patterns and spatiotemporal dynamics in chemical systems, Adv. Chem. Phys. 109 (2007) 435–513. [3] B. Rudovics, E. Barillot, P. Davies, E. Dulos, J. Boissonade, P. De Kepper, Experimental studies and quantitative modeling of Turing patterns in the (chlorine dioxide, iodine, malonic acid) reaction, J. Phys. Chem. A 103 (12) (1999) 1790–1800. [4] B. Rudovics, E. Dulos, P. De Kepper, Standard and nonstandard Turing patterns and waves in the CIMA reaction, Phys. Scr. 1996 (T67) (1996) 43. [5] L. Yang, I.R. Epstein, Oscillatory Turing patterns in reaction-diffusion systems with two coupled layers, Phys. Rev. Lett. 90 (17) (2003) 178303. [6] A. Coillet, I. Balakireva, R. Henriet, K. Saleh, L. Larger, J.M. Dudley, C.R. Menyuk, Y.K. Chembo, Azimuthal Turing patterns, bright and dark cavity solitons in Kerr combs generated with whispering-gallery-mode resonators, IEEE Photonics J. 5 (4) (2013) 6100409. [7] S. Kondo, T. Miura, Reaction-diffusion model as a framework for understanding biological pattern formation, Science 329 (5999) (2010) 1616–1620.

JID:YJCPH AID:109043 /FLA

18

[m3G; v1.261; Prn:4/11/2019; 13:14] P.18 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

[8] J.H. Cartwright, Labyrinthine Turing pattern formation in the cerebral cortex, J. Theor. Biol. 217 (1) (2002) 97–103. [9] M.D. Morales-Hernández, I.E. Medina-Ramírez, F. Avelar-González, J.E. Macías-Díaz, An efficient recursive algorithm in the computational simulation of the bounded growth of biological films, Int. J. Comput. Methods 9 (04) (2012) 1250050. [10] B. Pena, C. Perez-Garcia, Stability of Turing patterns in the Brusselator model, Phys. Rev. E 64 (5) (2001) 056213. [11] T. Biancalani, D. Fanelli, F. Di Patti, Stochastic Turing patterns in the Brusselator model, Phys. Rev. E 81 (4) (2010) 046215. [12] X. Tang, Y. Song, Bifurcation analysis and Turing instability in a diffusive predator-prey model with herd behavior and hyperbolic mortality, Chaos Solitons Fractals 81 (2015) 303–314. [13] T. Zhang, Y. Xing, H. Zang, M. Han, Spatio-temporal dynamics of a reaction-diffusion system for a predator–prey model with hyperbolic mortality, Nonlinear Dyn. 78 (1) (2014) 265–277. [14] F. Lutscher, A. Stevens, et al., Emerging patterns in a hyperbolic model for locally interacting cell systems, J. Nonlinear Sci. 12 (6) (2002) 619–640. [15] O.B. Isaeva, A.S. Kuznetsov, S.P. Kuznetsov, Hyperbolic chaos of standing wave patterns generated parametrically by a modulated pump source, Phys. Rev. E 87 (4) (2013) 040901. [16] E. Barbera, G. Consolo, G. Valenti, Spread of infectious diseases in a hyperbolic reaction-diffusion susceptible-infected-removed model, Phys. Rev. E 88 (5) (2013) 052719. [17] U.-I. Cho, B.C. Eu, Hyperbolic reaction-diffusion equations and chemical oscillations in the Brusselator, Phys. D, Nonlinear Phenom. 68 (3–4) (1993) 351–363. [18] M. Al-Ghoul, B.C. Eu, Hyperbolic reaction-diffusion equations, patterns, and phase speeds for the Brusselator, J. Phys. Chem. 100 (49) (1996) 18900–18910. [19] R. Eftimie, Hyperbolic and kinetic models for self-organized biological aggregations and movement: a brief review, J. Math. Biol. 65 (1) (2012) 35–75. [20] M. Wolfrum, The Turing bifurcation in network systems: collective patterns and single differentiated nodes, Phys. D, Nonlinear Phenom. 241 (16) (2012) 1351–1357. [21] J. Xu, G. Yang, H. Xi, J. Su, Pattern dynamics of a predator–prey reaction–diffusion model with spatiotemporal delay, Nonlinear Dyn. 81 (4) (2015) 2155–2163. [22] G. Consolo, C. Currò, G. Valenti, Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments, Appl. Math. Model. 43 (2017) 372–392. [23] V.E. Tarasov, Continuous limit of discrete systems with long-range interaction, J. Phys. A, Math. Gen. 39 (48) (2006) 14895. [24] V.E. Tarasov, G.M. Zaslavsky, Conservation laws and Hamilton’s equations for systems with long-range interaction and memory, Commun. Nonlinear Sci. Numer. Simul. 13 (9) (2008) 1860–1878. [25] R. Koeller, Applications of fractional calculus to the theory of viscoelasticity, J. Appl. Mech. (ISSN 0021-8936) 51 (1984) 299–307. [26] Y. Povstenko, Theory of thermoelasticity based on the space-time-fractional heat conduction equation, Phys. Scr. 2009 (T136) (2009) 014017. [27] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance, Phys. A, Stat. Mech. Appl. 284 (1) (2000) 376–384. [28] W.G. Glöckle, T.F. Nonnenmacher, A fractional calculus approach to self-similar protein dynamics, Biophys. J. 68 (1) (1995) 46–53. [29] V. Namias, The fractional order Fourier transform and its application to quantum mechanics, IMA J. Appl. Math. 25 (3) (1980) 241–265. [30] N. Su, P.N. Nelson, S. Connor, The distributed-order fractional diffusion-wave equation of groundwater flow: theory and application to pumping and slug tests, J. Hydrol. 529 (Part 3) (2015) 1262–1273. [31] V.G. Pimenov, A.S. Hendy, R.H. De Staelen, On a class of non-linear delay distributed order fractional diffusion equations, J. Comput. Appl. Math. 318 (2017) 433–443. [32] J.E. Macías-Díaz, Sufficient conditions for the preservation of the boundedness in a numerical method for a physical model with transport memory and nonlinear damping, Comput. Phys. Commun. 182 (12) (2011) 2471–2478. [33] J.E. Macías-Díaz, An explicit dissipation-preserving method for Riesz space-fractional nonlinear wave equations in multiple dimensions, Commun. Nonlinear Sci. Numer. Simul. 59 (2018) 67–87. [34] J.E. Macías-Díaz, On the solution of a Riesz space-fractional nonlinear wave equation through an efficient and energy-invariant scheme, Int. J. Comput. Math. 96 (2) (2018) 337–361. [35] A.A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys. 280 (2015) 424–438. [36] A.H. Bhrawy, M.A. Abdelkawy, A fully spectral collocation approximation for multi-dimensional fractional Schrödinger equations, J. Comput. Phys. 294 (2015) 462–483. [37] A. El-Ajou, O.A. Arqub, S. Momani, Approximate analytical solution of the nonlinear fractional KdV–Burgers equation: a new iterative algorithm, J. Comput. Phys. 293 (2015) 81–95. [38] F. Liu, P. Zhuang, I. Turner, V. Anh, K. Burrage, A semi-alternating direction method for a 2-D fractional FitzHugh–Nagumo monodomain model on an approximate irregular domain, J. Comput. Phys. 293 (2015) 252–263. [39] H. Ye, F. Liu, V. Anh, Compact difference scheme for distributed-order time-fractional diffusion-wave equation on bounded domains, J. Comput. Phys. 298 (2015) 652–660. [40] J.E. Macías-Díaz, A structure-preserving method for a class of nonlinear dissipative wave equations with Riesz space-fractional derivatives, J. Comput. Phys. 351 (2017) 40–58. [41] T. Langlands, B. Henry, S. Wearne, Turing pattern formation with fractional diffusion and fractional reactions, J. Phys. Condens. Matter 19 (6) (2007) 065115. [42] V. Gafiychuk, B. Datsko, Spatiotemporal pattern formation in fractional reaction-diffusion systems with indices of different order, Phys. Rev. E 77 (6) (2008) 066210. [43] B. Datsko, Y. Luchko, V. Gafiychuk, Pattern formation in fractional reaction–diffusion systems with multiple homogeneous states, Int. J. Bifurc. Chaos 22 (04) (2012) 1250087. [44] A. Mvogo, J.E. Macías-Díaz, T.C. Kofané, Diffusive instabilities in a hyperbolic activator-inhibitor system with superdiffusion, Phys. Rev. E 97 (3) (2018) 032129. [45] B.I. Henry, S.L. Wearne, Existence of Turing instabilities in a two-species fractional reaction-diffusion system, SIAM J. Appl. Math. 62 (3) (2002) 870–887. [46] Y. Nec, A. Nepomnyashchy, Turing instability in sub-diffusive reaction–diffusion systems, J. Phys. A, Math. Theor. 40 (49) (2007) 14687. [47] D. Jeong, Y. Choi, J. Kim, Modeling and simulation of the hexagonal pattern formation of honeycombs by the immersed boundary method, Commun. Nonlinear Sci. Numer. Simul. 62 (2018) 61–77. [48] D. Lacitignola, B. Bozzini, M. Frittelli, I. Sgura, Turing pattern formation on the sphere for a morphochemical reaction-diffusion model for electrodeposition, Commun. Nonlinear Sci. Numer. Simul. 48 (2017) 484–508. [49] X. Wang, W. Wang, G. Zhang, Vegetation pattern formation of a water-biomass model, Commun. Nonlinear Sci. Numer. Simul. 42 (2017) 571–584. [50] D. Prakasha, P. Veeresha, H.M. Baskonus, Two novel computational techniques for fractional Gardner and Cahn-Hilliard equations, Comput. Math. Methods 1 (2) (2019) e1021. [51] A.Q.M. Khaliq, X. Liang, K.M. Furati, A fourth-order implicit-explicit scheme for the space fractional nonlinear Schrödinger equations, Numer. Algorithms 75 (1) (2017) 147–172. [52] X. Liang, A.Q.M. Khaliq, H. Bhatt, K.M. Furati, The locally extrapolated exponential splitting scheme for multi-dimensional nonlinear space-fractional Schrödinger equations, Numer. Algorithms 76 (4) (2017) 939–958.

JID:YJCPH AID:109043 /FLA

[m3G; v1.261; Prn:4/11/2019; 13:14] P.19 (1-19)

J.E. Macías-Díaz / Journal of Computational Physics ••• (••••) ••••••

19

[53] K.M. Furati, M. Yousuf, A.Q.M. Khaliq, Fourth-order methods for space fractional reaction–diffusion equations with non-smooth data, Int. J. Comput. Math. 95 (6–7) (2018) 1240–1256. [54] Q.-J. Meng, D. Ding, Q. Sheng, Preconditioned iterative methods for fractional diffusion models in finance, Numer. Methods Partial Differ. Equ. 31 (5) (2015) 1382–1395. [55] J.E. Macías-Díaz, A.S. Hendy, R.H. De Staelen, A pseudo energy-invariant method for relativistic wave equations with Riesz space-fractional derivatives, Comput. Phys. Commun. 224 (2018) 98–107. [56] H. Ramos, Development of a new Runge-Kutta method and its economical implementation, Comput. Math. Methods 1 (2) (2019) e1016. [57] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, vol. 198, Elsevier, 1998. [58] J.E. Macías-Díaz, Numerical study of the transmission of energy in discrete arrays of sine-Gordon equations in two space dimensions, Phys. Rev. E 77 (1) (2008) 016602. [59] J.E. Macías-Díaz, A. Puri, On the transmission of binary bits in discrete Josephson-junction arrays, Phys. Lett. A 372 (30) (2008) 5004–5010. [60] W. Wang, Q.-X. Liu, Z. Jin, Spatiotemporal complexity of a ratio-dependent predator-prey system, Phys. Rev. E 75 (5) (2007) 051913. [61] M.D. Ortigueira, Riesz potential operators and inverses via fractional centred derivatives, Int. J. Math. Math. Sci. (2006). [62] C. Çelik, M. Duman, Crank–Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys. 231 (4) (2012) 1743–1750. [63] J.E. Macías-Díaz, Numerical simulation of the nonlinear dynamics of harmonically driven Riesz-fractional extensions of the Fermi–Pasta–Ulam chains, Commun. Nonlinear Sci. Numer. Simul. 55 (2018) 248–264. [64] K. Pen-Yu, Numerical methods for incompressible viscous flow, Sci. Sin. 20 (1977) 287–304.