Next Article in Journal
Carbon Fibers and Graphite as Pore-Forming Agents for the Obtention of Porous Alumina: Correlating Physical and Fractal Characteristics
Next Article in Special Issue
Zener Model with General Fractional Calculus: Thermodynamical Restrictions
Previous Article in Journal
Multifractal Characteristics of China’s Stock Market and Slump’s Fractal Prediction
Previous Article in Special Issue
A Brief Study on Julia Sets in the Dynamics of Entire Transcendental Function Using Mann Iterative Scheme
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Efficient Dissipation-Preserving Numerical Scheme to Solve a Caputo–Riesz Time-Space-Fractional Nonlinear Wave Equation

by
Jorge E. Macías-Díaz
1,2 and
Tassos Bountis
3,4,*
1
Department of Mathematics and Didactics of Mathematics, School of Digital Technologies, Tallinn University, Narva Rd. 25, 10120 Tallinn, Estonia
2
Departamento de Matemáticas y Física, Universidad Autónoma de Aguascalientes, Avenida Universidad 940, Ciudad Universitaria, Aguascalientes 20100, Mexico
3
Department of Mathematics, University of Patras, 26500 Patras, Greece
4
Center for Integrable Systems, P. G. Demidov Yaroslavl State University, 150003 Yaroslavl, Russia
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(9), 500; https://doi.org/10.3390/fractalfract6090500
Submission received: 5 August 2022 / Revised: 25 August 2022 / Accepted: 27 August 2022 / Published: 6 September 2022
(This article belongs to the Special Issue Feature Papers in Fractal and Fractional 2022–2023)

Abstract

:
For the first time, a new dissipation-preserving scheme is proposed and analyzed to solve a Caputo–Riesz time-space-fractional multidimensional nonlinear wave equation with generalized potential. We consider initial conditions and impose homogeneous Dirichlet data on the boundary of a bounded hyper cube. We introduce an energy-type functional and prove that the new mathematical model obeys a conservation law. Motivated by these facts, we propose a finite-difference scheme to approximate the solutions of the continuous model. A discrete form of the continuous energy is proposed and the discrete operator is shown to satisfy a conservation law, in agreement with its continuous counterpart. We employ a fixed-point theorem to establish theoretically the existence of solutions and study analytically the numerical properties of consistency, stability and convergence. We carry out a number of numerical simulations to verify the validity of our theoretical results.

1. Introduction

The investigation of mathematical models with fractional derivatives has become a fruitful area of research in recent decades. In fact, fractional models have been investigated in a wide range of scientific and engineering applications. For example, in the epidemiological sciences, such models have been employed to estimate the spreading of computer viruses [1], the propagation of measles in human populations [2], the transmission dynamics of varicella zoster virus [3], and even the propagation of the COVID-19 epidemic, including memory effects [4,5]. Other interesting applications of fractional calculus have focused on the design of ABC-fractional masks in image processing [6], the fractional modeling of blood ethanol concentration systems using real data [7], the investigation of three-species systems in biology using the Atangana–Baleanu fractional derivative [8], the mathematical modeling of the Chinese economic growth [9] and the description of cancer treatment, using a class of nonlinear fractional optimal control problems [10].
It is important to point out that the theory of fractional calculus itself has benefited tremendously from these contributions [11,12]. On the one hand, various new fractional operators have been proposed and analyzed in the literature [13,14], and new analytical results have been derived in each of these cases [15,16]. However, the broader significance of these new fractional operators in real-life applications is still an open issue. In fact, so far, only Riesz fractional operators have been adequately justified for their use in physical problems [17]. Indeed, fractional derivatives of the Riesz type have been proposed as a valid approach to the continuum limit of systems of particles with long-range interactions [18]. These operators have been studied mathematically, and various of their properties have been derived in the literature, including some of their associated variational structures [19]. These properties have already been employed in the investigation of thermodynamical systems, including the study of long-range interactions on the dynamics and statistics of 1D Hamiltonian lattices with on-site potentials [20] and supratransmission [21,22].
Classically, Caputo-like fractional operators have found potential applications in studying systems with memory effects. For instance, they have been used to model social and economic cycles of extreme events [23], elasticity in economic processes with memory [24] and the modeling of national economies in state-space [25]. It is worth pointing out that there are many types of fractional derivatives reported and analyzed in the specialized literature. However, only the Riesz and Caputo fractional operators seem to have physically relevant applications in several areas of science and engineering [26,27,28]. In particular, these operators have appeared in the extension of well-known parabolic or hyperbolic equations of mathematical physics, like fractional forms of the sine-Gordon equation [29,30], the nonlinear Klein–Gordon equation [31], the nonlinear Schrödinger equation [32], and the Klein–Gordon–Zakharov system [33]. It is worth mentioning that all of these models obey conservation laws which resemble a form of the mass or energy conservation of the equivalent systems in the integer-order case.
In the present work, we consider a fractional extension of the nonlinear Klein–Gordon equation with damping, which involves Caputo temporal derivatives and Riesz spatial derivatives. We show that this model satisfies a conservation law in differential form, and propose a discretization which satisfies the same law in the discrete domain. Our discretization leads to a finite-difference model which employs the L 1 scheme to approximate the Caputo temporal derivatives [34], and fractional-order centered differences to estimate the Riesz fractional derivatives in space [35]. Physically, this system is the continuum limit of a model describing particles with long-range interactions, memory effects and on-site potentials. For this reason, a theoretical study of this model and its discretization are of utmost importance. Various results from the standard literature will be used to that effect, including suitable forms of Gronwall’s inequality, and various estimates and properties for Caputo and Riesz fractional derivatives.
We must point out that there are various papers in the literature dealing with fractional hyperbolic systems which preserve the dissipation or conservation of energy property of classical systems of odes or pdes. Indeed, there are reports on fast conservative numerical algorithms for the coupled fractional Klein–Gordon–Schrödinger equation [36], conservative Fourier spectral methods of space-fractional Klein–Gordon–Schrödinger equations [37], energy-preserving finite-difference schemes for nonlinear fractional Schrödinger equations [38], mass-conservative Fourier spectral methods for solving the fractional nonlinear Schrödinger equation [39], an energy-preserving and efficient scheme for a double-fractional conservative Klein–Gordon–Zakharov system [40] and a conservative finite-difference scheme to solve a Riesz space-fractional Gross–Pitaevskii model [41]. However, in all of these works fractional derivatives (of the Riesz type) appear only in the spatial variables. It appears that the combined Caputo–Riesz time-space-fractional treatment of these systems involving conservation laws has not yet been investigated.
The present manuscript is one of the first works in the literature aiming to fulfill the tasks mentioned in the previous paragraph. More precisely, this is the first article in which a hyperbolic partial differential equation with fractional derivatives in both space and time is fully analyzed for a conserved quantity. The numerical model proposed in this work is capable of reflecting this property in the discrete domain as both theoretical and computational results will show. Moreover, the literature lacks reports in which the stability and convergence of finite-difference schemes for hyperbolic systems are theoretically derived. In that sense, this manuscript is also one of the first reports in the literature in which those studies are presented rigorously. It is important to mention that the numerical simulations confirm the validity of the mathematical results.

2. Preliminaries

Throughout this work, we suppose that p N represents the number of spatial dimensions which, for practical purposes, we may think of as p { 1 , 2 , 3 } , while Ω is a bounded, open, and connected domain of R p . More concretely, Ω is an open hyper-box of the form i = 1 p ( a i , b i ) , where a i , b i R and a i < b i . The real number T > 0 will denote a finite period of time, while Ω T = Ω × ( 0 , T ) . As usual, we let Ω and Ω T be the boundaries of Ω and Ω T , respectively. Throughout this work, α and β are real numbers in ( 1 , 2 ) , and define β 0 = β 1 . Moreover, all real functions with domain Ω T will be extended to the entire Euclidean domain R p × ( 0 , T ) by letting them be equal to zero in the exterior of the set Ω R p .
Definition 1
(Podlubny [42]).Let f : [ 0 , T ] R be a function, and suppose that β > 0 . Assume that m N { 0 } is a number with the property that m 1 < β m , and let t [ 0 , T ] . Whenever it exists, the Caputo fractional total derivative of f of order β at the time t is defined as
C D β f ( t ) = 1 Γ ( m β ) 0 t ( t s ) m 1 β f ( m ) ( s ) d s , if β Z , f ( m ) , if β = m Z .
Here, the symbol Γ represents the standard Gamma function. In turn, if  α > 0 satisfies α N , then the left and right Riesz fractional total derivatives of the function f of order α at the point x 0 are respectively defined as
d α f ( x ) d x α = 1 Γ ( m α ) d m d x m x ( x ξ ) m 1 α f ( ξ ) d ξ ,
d α f ( x ) d ( x ) α = ( 1 ) m Γ ( m α ) d m d x m x ( x ξ ) m 1 α f ( ξ ) d ξ .
If the left and right Riesz fractional total derivatives of f of order α at x 0 exist, then we define the respective Riesz total fractional derivatives of f by means of
d α f ( x ) d | x | α = 1 2 cos ( α π / 2 ) d α f ( x ) d x α + d α f ( x ) d ( x ) α .
In the case that α = m N , note that the Riesz fractional total derivative of f of order α is equal to the classical integer-order derivative of f of order m at x.
For the remainder, if n N , then we agree that I n = { 1 , , n } and I ¯ n = I n { 0 } . It is worth recalling that the Caputo fractional operator C D β satisfies various interesting properties. To start with, this operator (like the Riesz operator) is linear. Moreover, the Caputo derivative of a constant is equal to zero, and the following composition property is satisfied, for each β > 0 and n N { 0 } :
C D β f ( n ) ( t ) = C D β + n f ( t ) .
It is worthwhile mentioning that these last two properties are not shared by the Riesz fractional derivative [42].
Definition 2
(Podlubny [42]).Suppose that ϕ : Ω ¯ T R is a sufficiently smooth function, and let β and m be as in Definition 1. If  t [ 0 , T ] , then the temporal Caputo fractional partial derivative of ϕ of order β at ( x , t ) Ω T is given by
β ϕ ( x , t ) t β = 1 Γ ( m β ) 0 t ( t s ) m 1 β m ϕ ( x , s ) t m d s , if β Z , m ϕ ( x , t ) t m , if β = m Z .
If α is not an integer, i I p and x = ( x 1 , , x p ) Ω , we define respectively the left and right spatial Riesz fractional partial derivatives of ϕ of order α with respect to x i at ( x , t ) as
α ϕ ( x ) x i α = 1 Γ ( m α ) m x i m x ( x i ξ ) m 1 α ϕ ( x 1 , , x i 1 , ξ , x i + 1 , , x p ) d ξ ,
α ϕ ( x ) ( x i ) α = ( 1 ) m Γ ( m α ) m x i m x ( x i ξ ) m 1 α ϕ ( x 1 , , x i 1 , ξ , x i + 1 , , x p ) d ξ .
If these two quantities exist in the set of real numbers, we define the spatial Riesz fractional partial derivative of ϕ of order α with respect to x i at the point ( x , t ) Ω T by
α ϕ ( x , t ) | x i | α = 1 2 cos ( α π / 2 ) α ϕ ( x , t ) x i α + α ϕ ( x , t ) ( x i ) α .
As expected, if  α = m N , then the spatial Riesz partial derivative (9) is equal to the integer-order partial derivative of ϕ of order m with respect to x i at the point ( x , t ) . If all of these fractional partial derivatives exist, then we define the Riesz fractional Laplacian and gradient of ϕ of order α, respectively, as 
Δ α ϕ ( x , t ) = i = 1 p α ϕ ( x , t ) | x i | α ,
α ϕ ( x , t ) = α ϕ ( x , t ) | x 1 | α , α ϕ ( x , t ) | x 2 | α , , α ϕ ( x , t ) | x p | α .
The following are some useful properties on fractional derivatives of the Caputo and Riesz types. These properties will be used throughout what follows.
Lemma 1
(Alikhanov [43]).If f : [ 0 , T ] R is an absolutely continuous function and α ( 0 , 1 ) , then the following functional inequality holds:
1 2 C D β [ f 2 ( t ) ] f ( t ) C D β f ( t ) , t [ 0 , T ] .
For the remainder of this manuscript, we will suppose that ψ , χ : Ω ¯ R are sufficiently smooth functions. We will let γ 0 , and assume that F : R R is a continuously differentiable and non-negative function. Under these circumstances, the mathematical model under consideration here is the Caputo–Riesz time-space-fractional initial-boundary-value problem with homogeneous Neumann boundary conditions
β ϕ ( x , t ) t β + γ ϕ ( x , t ) t Δ α ϕ ( x , t ) + F ( ϕ ( x , t ) ) = 0 , ( x , t ) Ω T , such that ϕ ( x , 0 ) = ψ ( x ) , x Ω ¯ , ϕ ( x , 0 ) t = χ ( x ) , x Ω ¯ , ϕ ( x , t ) = 0 , ( x , t ) Ω × [ 0 , T ] .
It is obvious that this equation is a generalization of various particular integer-order models to the fractional scenario. Indeed, with  α = β = 2 , the resulting model becomes a hyperbolic system with damping. Under these circumstances, it reduces to the well-known sine-Gordon equation [44], if  F ( ϕ ) = 1 cos ϕ it leads to the linear Klein–Gordon equation of relativistic quantum mechanics [45], if  F ( ϕ ) = 1 2 ϕ 2 we get the nonlinear Klein–Gordon equation [46] for F ( ψ ) = 1 2 ! ϕ 2 1 4 ! ϕ 4 , and the double sine-Gordon equation [47] for F ( ϕ ) = 1 cos ϕ 1 2 cos ( 2 ϕ ) . In both the integer-order and the fully fractional cases, the fractional model (13) can be expressed equivalently as the coupled system of partial differential equations with initial-boundary conditions
β 0 υ ( x , t ) t β 0 + γ υ ( x , t ) Δ α ϕ ( x , t ) + F ( ϕ ( x , t ) ) = 0 , ( x , t ) Ω T , ϕ ( x , t ) t υ ( x , t ) = 0 , ( x , t ) Ω T , such that ϕ ( x , 0 ) = ψ ( x ) , x Ω ¯ , υ ( x , 0 ) = χ ( x ) , x Ω ¯ , ϕ ( x , t ) = 0 , ( x , t ) Ω × [ 0 , T ] .
It is important to point out that mathematical model (13) is obtained as the continuum-limit approximation for systems of particles with long-range interactions and memory effects. Indeed, it has been established rigorously that some discrete systems with long-range effects lead to systems of fractional differential equations when the limit between particles tends to zero. In those situations, the resulting fractional derivatives are Riesz operators in space [18,19]. Moreover, elasticity effects can be used to model systems with memory. Mathematically, this modeling can be achieved by means of Caputo-type fractional derivatives in the temporal variable [24].
Definition 3.
Define the set L x , p ( Ω ¯ T ) = { ϕ : Ω ¯ T R : ϕ ( · , t ) L p ( Ω ¯ ) , for each t [ 0 , T ] } and p [ 1 , ] . In the case that p [ 1 , ) and ϕ L x , p ( Ω ¯ T ) , we set
ϕ x , p = Ω ¯ | ϕ ( x , t ) | p d x 1 / p , t [ 0 , T ] .
If p = , then we let ϕ x , = inf { C 0 : | ϕ ( x , t ) | C for almost all x B ¯ } , for each t [ 0 , T ] . Obviously, the term ϕ x , p is a function of t [ 0 , T ] in any case. Moreover, for each pair ϕ , υ L x , 2 ( Ω ¯ ) , define the following function of t:
ϕ , υ x = Ω ¯ ϕ ( x , t ) υ ( x , t ) d x , t [ 0 , T ] .
Lemma 2.
Let ϕ , υ : Ω ¯ T R be sufficiently smooth functions, and  α ( 1 , 2 ) . If  i I p , then
α ϕ | x i | α , υ x = ϕ , α υ | x i | α x = α / 2 ϕ | x i | α , α / 2 υ | x i | α / 2 x .
Theorem 1
(Conservation law).If ϕ and υ are solutions of (14), then E ( t ) = 0 , where
E ( t ) = β 0 υ t β 0 , υ x + γ υ x , 2 2 + 1 2 d d t α / 2 ϕ x , 2 2 + d d t F ( ϕ ) , 1 x , t [ 0 , T ] .
Proof. 
We will provide only a sketch of the proof, since most of the details result from algebraic simplifications and taking derivatives. Beforehand, it is easy to check that the following identities are satisfied:
Δ α ϕ , ϕ t = i = 1 p α ϕ | x i | α , ϕ t = i = 1 p α / 2 ϕ | x i | α / 2 , t α / 2 ϕ | x i | α / 2 = 1 2 i = 1 p t α / 2 ϕ | x i | α / 2 2 = 1 2 d d t α / 2 ϕ x , 2 2 ,
and
F ( ϕ ) , ϕ t = Ω ¯ F ( ϕ ) ϕ t d x = Ω ¯ F ( ϕ ) t d x = d d t F ( ϕ ) , 1 x .
Consider the first equation of (14), take the inner product on both sides with respect to υ , and recall that the second equation of that system assures that υ is actually the first-order partial derivative of ϕ with respect to t. Finally, use the above identities and recall that β = β 0 + 1 to reach the statement of the Theorem.    □
As an easy consequence of our last result and a direct integration over the interval [ 0 , t ] [ 0 , T ] , we readily reach the following conservation law, valid for each t [ 0 , T ] :
0 t β 0 υ t β 0 , υ x d s + 0 t γ υ x , 2 2 d s + 1 2 α / 2 ϕ ( · , t ) x , 2 2 + F ( ϕ ( · , t ) ) , 1 x = 1 2 α / 2 ψ x , 2 2 + F ( ψ ) , 1 x .
Corollary 1.
If ϕ and υ are solutions of (14), then
1 2 β 0 υ x , 2 2 t β 0 + 1 2 d d t α / 2 ϕ x , 2 2 + d d t F ( ϕ ) , 1 x 0 , t [ 0 , T ] .
Proof. 
The result is a direct consequence from Lemma 1 and Theorem 1.    □
This corollary and the conservation law (21) yield that the following property is satisfied for each t [ 0 , T ] :
1 2 0 t β 0 υ x , 2 2 t β 0 d s + 1 2 α / 2 ϕ ( · , t ) x , 2 2 + F ( ϕ ( · , t ) ) , 1 x 1 2 α / 2 ψ x , 2 2 + F ( ψ ) , 1 x .

3. Numerical Model

At this stage, we introduce the discrete systems needed to approximate the solutions of the mathematical model (14), introducing suitable notation to describe the numerical model and its discrete functionals. Our approach will be based on the use of finite differences. In this way, we will suppose that h i > 0 and τ > 0 , for each i I p . Additionally, for each i I p , assume also that N = T / τ and M i = ( b i a i ) / h i are both positive integers.
Let us take regular partitions of the spatial intervals [ a i , b i ] and the temporal interval [ 0 , T ] , given, respectively, by 
x i , j i = a i + j i h i , i I p , j i I ¯ M i ,
t n = n τ , n I ¯ N ,
and define t n + 1 2 = ( n + 1 2 ) τ , for each n I ¯ N 1 . As is customary in our investigation [48], we introduce the collections J = i = 1 p I M i 1 and J ¯ = i = 1 p I ¯ M i . Additionally, the notation J will represent the boundary of the grid J ¯ , which is the set of all multi-indexes j J ¯ with the property that x j Ω . Here, we agree tacitly that x j = ( x 1 , j 1 , , x p , j p ) , for each j = ( j 1 , , j p ) J ¯ .
The symbols Φ j n and Υ j n represent numerical approximations to the exact value of ϕ j n = ϕ ( x j , t n ) and υ j n = υ ( x j , t n ) , respectively, for each j J ¯ and n I ¯ N . For simplicity, we set ϕ n = ( ϕ j n ) j J ¯ and υ n = ( υ j n ) j J ¯ . In a similar fashion, we let Φ n = ( Φ j n ) j J ¯ and Υ n = ( Υ j n ) j J ¯ , and agree that Φ = ( Φ n ) n I ¯ N and Υ = ( Υ n ) n I ¯ N .
Definition 4.
Let Ψ represent any one of the symbols Φ or Υ, and let j J and n I N . We define the discrete linear operators
δ τ Ψ j n = Ψ j n + 1 Ψ j n τ ,
μ τ Ψ j n = Ψ j n + 1 + Ψ j n 2 .
Moreover, if  F : R R is a continuously differentiable function, j J and n I N , then we let
δ Ψ n F ( Ψ j n ) = F ( Ψ j n + 1 ) F ( Ψ j n ) Ψ j n + 1 Ψ j n , if Ψ j n + 1 Ψ j n , F ( Ψ j n ) , o t h e r w i s e .
It is worth pointing out that there are no general approaches to deal with the nonlinear term in (13). In our present work, we employ the discretization formula (28). This disretization of the nonlinear term has been chosen because it is useful to preserve variational properties of the solutions for the continuous system. Ultimately, these variational arguments are used in the form of the energy method, in order to establish the stability and convergence for the numerical scheme.
Definition 5
(Zhuang and Liu [49]).Let τ > 0 , assume that f : [ 0 , ) R is sufficiently regular, and let β ( 0 , 1 ) . Let us define t k = k τ , for each k N { 0 } , and define the coefficients
b m ( β ) = τ 1 β Γ ( 2 β ) ( m + 1 ) 1 β m 1 β , m N { 0 } .
Recall that the L 1 scheme of the function f of order β at the time t k is given by
δ τ ( β ) f ( t k + 1 ) = 1 τ m = 0 k b k m ( β ) f ( t m + 1 ) f ( t m ) , k N .
Various remarks are in order. To start with, notice that the L 1 scheme is not defined when β = 1 since the Gamma function is not defined when the argument is equal to zero. On the other hand, it is well known [50] that, if  f C 2 ( [ 0 , ) ] has all its derivatives up to order 2 bounded on [ 0 , ) , then the following consistency property is satisfied:
δ τ ( β ) f ( t k ) = d β f ( t k ) d t β + O ( τ 2 β ) , k N .
Moreover, if  β ( 1 , 2 ) , then we define the discrete linear operator δ τ ( β ) = δ τ ( β 0 ) δ τ . Obviously, ∘ in this case represents the operation of composition. Under these circumstances, δ τ ( β ) provides a consistent approximation for the Caputo fractional derivative operator of order β , for functions which belong to C 3 ( [ 0 , ) ] and possess derivatives up to order 3 that are bounded on [ 0 , ) . In such a case, it is possible to check that the consistency order is τ 3 β (see [50]).
Definition 6
(Ortigueira [51]).Assume that α > 1 , and define the constants
g k ( α ) = ( 1 ) k Γ ( α + 1 ) Γ ( α 2 k + 1 ) Γ ( α 2 + k + 1 ) , k Z .
In the case that α ( 1 , 2 ] , it is easy to check that g 0 ( α ) 0 , and that g k ( α ) = g k ( α ) < 0 , for each k N . Moreover, using results from [52], it follows that
k = g k ( α ) = 0 .
Definition 7
(Ortigueira [51]).Suppose that h > 0 and that f : R R is any function. Then the fractional centered difference of f of order α of f at x R is defined as
Δ h α f ( x ) = k = g k ( α ) f ( x k h ) .
If α ( 1 , 2 ] and f C 5 ( R ) has all of its derivatives up to order five in L 1 ( R ) , then [52]
1 h α Δ h α f ( x ) = α f ( x ) | x | α + O ( h 2 ) , x R .
Definition 8.
Let i I p , j J ¯ , α ( 1 , 2 ] and β ( 0 , 1 ) . We define the discrete operators
δ τ ( β ) Φ j n + 1 = m = 0 n b n m ( β ) Φ ( x j , t m + 1 ) Φ ( x j , t m ) τ ,
δ x i ( α ) Φ j n = 1 h i α k = 0 M i g j i k ( α ) Φ ( x 1 , j 1 , , x i 1 , j i 1 , x i , k , x i + 1 , j i + 1 , , x p , j p , t n ) .
In view of the preceding discussion, the first operator provides an approximation to the Caputo partial derivative of ϕ with respect to t at the point ( x j , t n + 1 ) , with a consistency order equal to 2 β . In turn, the second operator yields a quadratically consistent approximation to the Riesz fractional partial derivatives of ϕ with respect to x i at the point x j and time t n . Moreover, notice that if β ( 1 , 2 ) , then we will write β τ ( β ) = δ τ ( β 0 ) δ τ . Obviously, in such case, one may readily check that
δ τ ( β ) Φ j n + 1 = m = 0 n b n m ( β ) Φ ( x j , t m + 1 ) 2 Φ ( x j , t m ) + Φ ( x j , t m 1 ) τ 2 = b 0 ( β ) Φ ( x j , t n + 1 ) 2 Φ ( x j , t n ) + Φ ( x j , t n 1 ) τ 2 + m = 0 n 1 b n m ( β ) Φ ( x j , t m + 1 ) 2 Φ ( x j , t m ) + Φ ( x j , t m 1 ) τ 2 ,
where
b m ( β ) = τ 2 β Γ ( 3 β ) ( m + 1 ) 2 β m 2 β , m N { 0 } .
Definition 9.
Using the same conventions as in the previous definition, we may approximate the spatial Laplacian and gradient operators, respectively, by 
Δ h ( α ) Φ j n = i = 1 p δ x i ( α ) Φ j n , ( j , n ) J ¯ × I ¯ N ,
h ( α ) Φ j n = δ x i ( α ) Φ j n , δ x i ( α ) Φ j n , , δ x i ( α ) Φ j n , ( j , n ) J ¯ × I ¯ N .
Adopting the nomenclature and the discrete operators introduced in the present section, the finite-difference method that approximates the solutions of the system (14) is given by the discrete initial-boundary-value problem
δ τ ( β 0 ) Υ j n + γ δ τ Φ j n μ τ Δ h ( α ) Φ j n + δ Φ n F ( Φ j n ) = 0 , δ τ Φ j n Υ j n = 0 , such that Φ j 0 = ψ ( x j ) , j J ¯ , Υ j n = χ ( x j ) , j J ¯ , Φ j n = Υ j n = 0 , ( j , n ) J × I ¯ N .
It is obvious that this system is a two-step, implicit algebraic model. The first equation is a nonlinear system in view of the presence of the last term at the left-hand side. Meanwhile, the second equation yields a linear system which can be described easily in vector form. The overall system is nonlinear, and requires the computational implementation of a suitable method to solve the above equations. In our simulations at the end of this work, we will employ a version of the Newton–Raphson technique to obtain its solution.
Note that the above recursive system of finite-difference equations can be rewritten equivalently as a single discrete iterative equation
δ τ ( β ) Φ j n + γ δ τ Φ j n μ τ Δ h ( α ) Φ j n + δ Φ n F ( Φ j n ) = 0 .
Here, we define δ τ ( β ) Φ j n = ( δ τ ( β 0 ) δ τ ) Φ j n , for each ( j , n ) J ¯ × I ¯ N 1 .
Definition 10.
For the remainder of this manuscript, we set h = ( h 1 , h 2 , , h p ) , and define the constant h * = h 1 h 2 h p . Fix now R h = { x j : j J ¯ } , and let V h be the collection of real-valued functions defined on R h which vanish on J . Let us introduce the functions · p : V h R and · , · : V h × V h R by
Φ p = h * j J ¯ | Φ j | p 1 / p , Φ V h ,
Φ , Ψ = h * j J ¯ Φ j Ψ j , Φ , Ψ V h .
Notice that the first function is a norm, while the second is an inner product. Moreover, let · : V h R denote the infinity norm in V h , which is given by Φ = max { | Φ j | : j J ¯ } , for each Φ V h . With this notation, it is easy to check that ϕ n , υ n , Φ n , Υ n V h , for each n I ¯ N .
The following is a discrete version of Alikhanov’s inequality in Lemma 1.
Lemma 3
(Alikhanov [53]).If Φ V h and β ( 0 , 1 ) , then 1 2 δ τ ( β ) Φ 2 2 δ τ ( β ) Φ , Φ .
Lemma 4
(Macías-Díaz [54]).Let α ( 1 , 2 ] , and suppose that i I p . If  Φ V h , then the inequality δ x i ( α ) Φ 2 2 4 ( g 0 ( α ) h * h i α ) 2 Φ 2 2 is satisfied. Moreover, if  Φ , Ψ V h , then
δ x i ( α ) Φ , Ψ = Φ , δ x i ( α ) Ψ = δ x i ( α / 2 ) Φ , δ x i ( α / 2 ) Ψ .
Additionally, h ( α ) Φ 2 g h ( α ) Φ 2 , where
g h ( α ) = 2 g 0 ( α ) h * p i = 1 p h i 2 α 1 / 2 .
As our first theoretical result, we would like to prove the solvability of the finite-difference method (42). To that end, we will need the following classical result from the literature.
Lemma 5
(Brouwer fixed-point theorem [55]).Assume that ( H , · , · ) is an inner-product real or complex vector space of finite dimension, and denote with · the induced Euclidean norm. Suppose that g : H H is a continuous function, and assume that there exists some β > 0 , such that Re g ( z ) , z > 0 , for each z H with z = β . Then there exists some z * H with the properties that g ( z * ) = 0 and z * β .
Let Φ = ( Φ n ) n I ¯ N and Υ = ( Υ n ) n I ¯ N be two finite sequences in V h , let Ψ V h and assume that F : R R is a differentiable function. For the sake of convenience, we introduce the discrete operator δ Ψ , Φ n F ( Φ n ) , which is defined for each ( j , n ) J ¯ × I ¯ N by
δ Ψ , Φ n F ( Φ j n ) = F ( Ψ j ) F ( Φ j n ) Ψ j Φ j n , if Ψ j Φ j n , F ( Φ j n ) , otherwise .
Notice that when the function F is differentiable and its derivative is bounded, then the mean value theorem guarantees that there exists a constant K 0 which is independent of Ψ , such that δ Ψ , Φ n F ( Φ n ) K , for each n I ¯ N .
Now, for each ( j , n ) J ¯ × I ¯ N 1 , we define the discrete operator
δ τ , Φ n + 1 ( β ) Ψ j = b 0 ( β ) Ψ j 2 Φ j n + Φ j n 1 τ 2 + m = 0 n 1 b n m ( β ) Φ j m + 1 2 Φ j m + Φ j m 1 τ 2 , Ψ V h .
After some easy calculations involving the Cauchy–Schwartz inequality, it is easy to check that there exists a non-negative real-valued function G n + 1 = G n + 1 ( Φ ) which only depends on Φ 0 , Φ 1 , , Φ n + 1 , with the property that the following inequality is satisfied:
δ τ , Φ n + 1 ( β ) Ψ , Ψ b 0 ( β ) τ 2 Ψ 2 Ψ 2 G n + 1 , n I ¯ N 1 , Ψ V h .
All these facts will be employed in the proof of the theorem that follows.
Theorem 2
(Existence of solutions).Let F : R R be a differentiable function with bounded derivative in R . Then the finite-difference method (42) is solvable for any set of initial conditions if the following inequality is satisfied:
2 τ β Γ ( 3 β ) + 2 γ τ > g h ( α ) .
Proof. 
We will prove this result using induction. Notice firstly that the initial approximations are defined explicitly by the initial data of (42). So, let us assume that the approximations Φ n and Υ n have already been calculated, for some n I ¯ N 1 . Let us define the function g : V h V h by
g ( Ψ ) = δ τ , Φ n ( β ) Ψ + γ τ ( Ψ Φ n ) 1 2 Δ h ( α ) ( Ψ + Φ n ) + δ Ψ , Φ n F ( Φ n ) , Ψ V h ,
Next, we use the Cauchy–Schwartz and Young’s inequalities, the remarks preceding the statement of this theorem and Lemma 4 together with some algebraic regrouping to derive the following inequalities, for each Ψ V h :
γ τ Ψ Φ n , Ψ γ τ Ψ 2 Ψ 2 Φ n 2 ,
1 2 Δ h ( α ) ( Ψ + Φ n ) , Ψ 1 2 Ψ 2 g h ( α ) Ψ 2 + Δ h ( α ) Φ n 2 ,
δ Ψ , Φ n F ( Φ n ) , Ψ K Ψ 2 .
Adding these inequalities and regrouping, we obtain that g ( Ψ ) , Ψ Ψ 2 ζ Ψ 2 η , where the constants ζ and η are given by
ζ = b 0 ( β ) τ 2 + γ τ 1 2 g h ( α ) = 2 b 0 ( β ) + 2 τ γ τ 2 g h ( α ) 2 τ 2 ,
η = b 0 ( β ) τ 2 G n + 1 + γ τ Φ n 2 2 + 1 2 Δ h ( α ) Φ n 2 2 + K ,
by virtue of the fact that the inequality (51) is satisfied by hypothesis, ζ > 0 . On the other hand, notice that η is a number which depends only on the model parameters and the approximations Φ 0 , Φ 1 , , Φ n . It follows that g ( Ψ ) , Ψ 0 whenever Ψ V h satisfies Ψ 2 = η / ζ . Brouwer’s fixed-point theorem assures now that there exists some Φ n + 1 V h such that g ( Φ n + 1 ) = 0 . Indeed, just take λ = η / ζ in Lemma 5. The conclusion of this theorem readily follows now by mathematical induction.    □
The following result is a discrete version of Theorem 1.
Theorem 3
(Discrete conservation law).If Φ and Υ are solutions of the discrete model (42), then E n = 0 for each n I ¯ N 1 , where
E n = δ τ ( β 0 ) Υ n , Υ n + γ δ τ Φ n 2 2 + 1 2 δ τ h ( α / 2 ) Φ n 2 2 + δ τ F ( Φ n ) , 1 .
Proof. 
To establish this result, it is important to notice first that the following identities are satisfied, for each n I ¯ N 1 :
μ τ Δ h ( α ) Φ n , δ τ Φ n = i = 1 p δ x i ( α ) μ τ Φ n , δ τ Φ n = i = 1 p μ τ δ x i ( α / 2 ) Φ n , δ τ δ x i ( α / 2 ) Φ n = 1 2 τ i = 1 p δ x i ( α / 2 ) Φ n + 1 2 2 δ x i ( α / 2 ) Φ n 2 2 = 1 2 i = 1 p δ τ δ x i ( α / 2 ) Φ n 2 2 = 1 2 δ τ h ( α / 2 ) Φ n 2 2 ,
and
δ Φ F ( Φ n ) , δ τ Φ n = 1 τ F ( Φ n + 1 ) F ( Φ n ) , 1 = δ τ F ( Φ n ) , 1 .
Take now the inner product of the first equation of the finite-difference scheme (42) with δ τ Φ n , and use the identities above. Finally, use the second equation of (42) to calculate the first term, and thus one readily reaches the identity (58).    □
Let Ψ j = ψ ( x j ) and X j = χ ( x j ) , for each j J ¯ . Multiplying both sides of (58) by τ , summing over the first n time steps and using telescoping sums, we readily obtain the following identity, which is the discrete version of (21), for each n I ¯ N 1 :
k = 0 n δ τ ( β 0 ) Υ k , Υ k τ + γ k = 0 n δ τ Φ k 2 2 τ + 1 2 h ( α / 2 ) Φ n + 1 2 2 + F ( Φ n + 1 ) , 1 = 1 2 h ( α / 2 ) Ψ 2 2 + F ( X ) , 1 .
The next result is a straightforward consequence of Lemma 3 and Theorem 3. It is worth pointing out that, additionally, it is a fully discrete form of Corollary 1.
Corollary 2.
If Φ and Υ satisfy (42), then
1 2 δ τ ( β 0 ) Υ n 2 2 + 1 2 δ τ h ( α / 2 ) Φ n 2 2 + δ τ F ( Φ n ) , 1 0 , n I ¯ N 1 .
As a consequence of this last corollary, we readily obtain the dissipation rule
1 2 τ k = 0 n δ τ ( β 0 ) Υ k 2 2 + 1 2 h ( α / 2 ) Φ n + 1 2 2 + F ( Φ n + 1 ) , 1 1 2 h ( α / 2 ) Ψ 2 2 + F ( X n + 1 ) , 1 , n I ¯ N 1 .

4. Numerical Properties

The purpose of the present section is to rigorously establish the most important properties of the finite-difference scheme (42). Concretely, we prove theoretically the consistency, stability, and convergence of the numerical model. It is worth pointing out that the uniqueness of the solutions of the discrete system (42) will follow as a consequence of the stability of the finite-difference scheme.
First, we prove the consistency of the above numerical model and, to that end, we need to introduce the following continuous differential operators:
L ( x , t ) = β 0 υ ( x , t ) t β 0 + γ ϕ ( x , t ) t Δ α ϕ ( x , t ) + F ( ϕ ( x , t ) ) , ( x , t ) Ω T ,
M ( x , t ) = ϕ ( x , t ) t υ ( x , t ) , ( x , t ) Ω T .
To simplify matters, we define L j n = L ( x j , t n + 1 2 ) and M j n = M ( x j , t n + 1 2 ) . Using this notation, we let L n = ( L j n ) j J ¯ and M n = ( M j n ) j J ¯ , for each n I ¯ N 1 . Moreover, we set L = ( L n ) n I ¯ N 1 and M = ( M n ) n I ¯ N 1 .
Let us introduce the discrete difference operators
L j n = δ τ ( β 0 ) υ j n + γ δ τ ϕ j n μ τ Δ h ( α ) ϕ j n + δ ϕ F ( ϕ j n ) , ( j , n ) J × I ¯ N 1 ,
M j n = δ τ ϕ j n υ j n , ( j , n ) J × I ¯ N 1 .
Now set L n = ( L j n ) j J ¯ and M n = ( M j n ) j J ¯ , for each i I ¯ N 1 , and let L = ( L n ) n I ¯ N 1 and M = ( M n ) n I ¯ N 1 .
Lemma 6
(Sun and Wu [56]).If g C 2 ( [ 0 , T ] and 1 < β < 2 , then
0 t n g ( t ) d t ( t n t ) β 1 1 τ b 0 ( β ) g ( t n ) k = 1 n 1 ( b n k 1 ( β ) b n k ( β ) ) g ( t k ) b n 1 ( β ) g ( t 0 ) 1 2 β 2 β 12 + 2 3 β 3 β ( 1 + 2 1 β ) max 0 t t n g ( t ) τ 3 β ,
where
b k ( β ) = τ 2 β 2 β ( k + 1 ) 2 β k 2 β , k N { 0 } .
The previous result will be essential below in proving the consistency of the finite-difference scheme (42). Moreover, we will require the fact that δ τ ( β ) = δ τ ( β 0 ) δ τ , for each 1 < β < 2 . Here, the symbol ∘ represents composition of operators.
Theorem 4
(Consistency).Suppose that ( ϕ , υ ) is a solution of (14), and that ϕ C x , t 5 , 3 ( Ω ¯ T ) and υ C t 2 ( Ω T ) . If  τ < 1 , then there exists a constant C 0 which is independent of h and τ, with the property that
max | L L | , | M M | C ( τ 3 β + h 2 2 ) .
Proof. 
As is customary, the proof of this property hinges on the use of the regularity assumptions on the solutions of the continuous problem (14), the smoothness of F, Taylor’s theorem, and the previous lemma together with the consistency properties of the L 1 -scheme and the fractional-order centered differences. As a first step, notice that the regularity of ϕ and υ assure that there exist constants C 1 , C 2 , C 3 i , C 4 , C 5 R + { 0 } independent of h and τ , such that the following estimates hold, for each i I p and ( j , n ) J × I ¯ N 1 :
β 0 υ ( x j , t n + 1 2 ) t β 0 μ τ δ τ ( β 0 ) υ j n C 1 τ 3 β ,
ϕ ( x j , t n + 1 2 ) t δ τ ϕ j n C 2 τ 2 ,
α ϕ ( x j , t n + 1 2 ) | x i | α μ τ δ x i ( α ) ϕ j n C 3 i ( τ 2 + h i 2 ) ,
F ( ϕ ( x j , t n + 1 2 ) δ ϕ F ( ϕ j n ) C 4 τ 2 ,
υ ( x j , t n + 1 2 ) μ τ υ j n C 5 τ 2 .
As a consequence of inequality () and the triangle inequality, there exists a constant C 3 0 which is independent of h and τ such that, for each ( j , n ) J × I ¯ N 1 , the following inequality is satisfied:
Δ α ϕ ( x j , t n + 1 2 ) μ τ Δ h ( α ) ϕ j n C 3 ( τ 2 + h 2 2 ) .
Recall that τ < 1 and that 3 β 2 . The conclusion of this result follows now by the triangle inequality and letting C = max { C 1 * , C 2 * } , where the nonnegative constants C 1 and C 2 are given by C 1 * = C 1 + γ C 2 + C 3 + C 4 and C 2 * = C 2 + C 5 , respectively. It is obvious that this construction guarantees that C is independent of h and τ , as required. Indeed, notice beforehand that the hypothesis τ < 1 is equivalent to τ β 1 < 1 . In turn, this fact implies that τ 2 < τ 3 β . Observe then that the triangle inequality yields
| L j n L j n | β 0 υ ( x j , t n + 1 2 ) t β 0 μ τ δ τ ( β 0 ) υ j n + γ ϕ ( x j , t n + 1 2 ) t δ τ ϕ j n + Δ α ϕ ( x j , t n + 1 2 ) μ τ Δ h ( α ) ϕ j n + F ( ϕ ( x j , t n + 1 2 ) δ ϕ F ( ϕ j n ) C 1 τ 3 β + γ C 2 τ 2 + C 3 ( τ 2 + h 2 2 ) + C 4 τ 2 C 1 * ( τ 3 β + h 2 2 ) ,
for each ( j , n ) J × I ¯ N 1 . Obviously, the constant C 2 * is also independent of h and τ , and is obtained in similar fashion.    □
Definition 11.
The Mittag-Leffler function with one parameter β ( 0 , 1 ) is the complex function given for every z C by
E β ( z ) = k = 0 z k Γ ( 1 + k β ) .
It is well known that the Mittag-Leffler function provides an extension of various well-known complex functions. For example, E 1 ( z ) = e z , E 2 ( z ) = cosh ( z ) and E 1 2 ( z ) = e z 2 erfc ( z ) , for each z C . Here, the function erfc is the error function complement [42].
Our next goal is to establish the property of stability of the finite-difference scheme (42). To that end, we will need the following results from the literature:
Lemma 7
(Hendy and Macías-Díaz [57]).Let 0 < β < 1 . Suppose that ( ω n ) n = 0 and ( g n ) n = 0 are non-negative sequences, and let C 1 and C 2 be non-negative constants. If  C 1 > 0 and
δ τ ( β ) ω n C 1 ω n + C 2 ω n 1 + g n , n I N ,
then there exists a positive constant τ * with the property that, for each 0 < τ τ * , the following inequality is satisfied, for each n I ¯ N :
ω n 2 ω 0 + t n β Γ ( 1 + β ) max 0 n 0 n g n 0 E β ( 2 C 0 t n β ) .
Here,
C 0 = C 1 + C 2 2 2 1 β .
Lemma 8
(Macías-Díaz [54]).Let F C 2 ( R ) and F L ( R ) , and suppose that ( Φ n ) n I ¯ N , ( Φ ˜ n ) n I ¯ N and ( R n ) n I ¯ N are finite sequences in V h . Let us assume that ϵ n = Φ n Φ ˜ n and
F ˜ Φ , Φ ˜ n = δ Φ n F ( Φ n ) δ Φ ˜ n F ( Φ ˜ n ) , n I ¯ N .
Then there exists a constant C 1 0 which depends only on F, with the property that
2 τ n = 1 k R n F ˜ Φ , Φ ˜ n , δ τ ϵ n 2 τ n = 0 k R n 2 2 + C 1 ϵ 0 2 2 + C 1 τ n = 0 k δ τ ϵ n 2 2 , k I ¯ N .
Theorem 5
(Stability).Let F C 2 ( R ) and F , F L ( R ) , and suppose that (51) holds. Assume that ( ψ , χ ) and ( ψ ˜ , χ ˜ ) are initial conditions for the discrete model (42), with solutions ( Φ , Υ ) and ( Φ ˜ , Υ ˜ ) , respectively. Let ϵ n = Φ n Φ ˜ n , for each n I ¯ N , and let C 1 be as in Lemma 8. Then there exists τ * > 0 such that, for  0 < τ τ * , the following is satisfied, for each n I ¯ N 1 :
ϵ n + 1 ϵ n 2 2 2 ϵ 1 ϵ 0 2 2 + T 1 + β Γ ( 1 + β ) ( h ( α / 2 ) ϵ 0 2 2 + C 1 ϵ 0 2 2 ) E β ( 2 C 1 T β ) .
Proof. 
Since inequality (51) is satisfied, then the solutions ( Φ , Υ ) and ( Φ ˜ , Υ ˜ ) corresponding to the initial conditions ( ψ , χ ) and ( ψ ˜ , χ ˜ ) , respectively, do exist. Moreover, the first solution satisfies the finite-difference scheme (42), while the second satisfies
δ τ ( β 0 ) Υ ˜ j n + γ δ τ Φ ˜ j n μ τ Δ h ( α ) Φ ˜ j n + δ Φ ˜ n F ( Φ ˜ j n ) = 0 , δ τ Φ ˜ j n Υ ˜ j n = 0 , such that Φ ˜ j 0 = ψ ˜ ( x j ) , j J ¯ , Υ ˜ j n = χ ˜ ( x j ) , j J ¯ , Φ ˜ j n = Υ ˜ j n = 0 , ( j , n ) J × I ¯ N .
Let us define the vectors ϵ n = Φ n Φ ˜ n and ζ n = Υ n Υ ˜ n , for each n I ¯ N . Subtracting equations (85) from those of (42), it is possible to check that the pair ( ϵ , ζ ) satisfies the following discrete initial-boundary-value problem
δ τ ( β 0 ) ζ j n + γ δ τ ϵ j n μ τ Δ h ( α ) ϵ j n + F ˜ Φ , Φ ˜ n = 0 , δ τ ϵ j n ζ j n = 0 , such that ϵ j 0 = ψ ( x j ) ψ ˜ ( x j ) , j J ¯ , ζ j n = χ ( x j ) χ ˜ ( x j ) , j J ¯ , ϵ j n = ζ j n = 0 , ( j , n ) J × I ¯ N ,
where F ˜ Φ , Φ ˜ n = δ Φ n F ( Φ n ) δ Φ ˜ n F ( Φ ˜ n ) . Take now the inner product of the first identity in (86) with δ τ ϵ n , substitute the second identity rearrange terms algebraically and apply Lemma 3, to obtain
1 2 δ τ ( β 0 ) δ τ ϵ n 2 2 + 1 2 δ τ h ( α / 2 ) ϵ n 2 2 1 2 δ τ ( β 0 ) δ τ ϵ n 2 2 + γ δ τ ϵ n 2 2 + 1 2 δ τ h ( α / 2 ) ϵ n 2 2 δ τ ( β 0 ) μ τ ζ n , μ τ ζ n + γ δ τ ϵ n 2 2 + 1 2 δ τ h ( α / 2 ) ϵ n 2 2 = F ˜ Φ , Φ ˜ n , δ τ ϵ n F ˜ Φ , Φ ˜ n , δ τ ϵ n , n I ¯ N .
Now, let k I N , and take the sum on both ends of this chain of identities and inequalities for all n from 0 to k 1 . Multiply then by 2 τ 2 , apply the formula for telescoping sums, use Lemma 8 with R n = 0 to bound from above and rearrange terms algebraically to reach
δ τ ( β 0 ) ω k 1 τ 2 n = 0 k 1 δ τ ( β 0 ) δ τ ϵ n 2 2 + τ h ( α / 2 ) ϵ k 2 2 τ h ( α / 2 ) ϵ 0 2 2 + 2 τ 2 n = 0 k 1 F ˜ Φ , Φ ˜ n , δ τ ϵ n g + C 1 ω k 1 .
Here, the constant C 1 is the one from Lemma 8. It depends only on the function F, and
ω k = n = 0 k ϵ n + 1 ϵ n 2 2 , k I ¯ N 1 ,
g = T ( h ( α / 2 ) ϵ 0 2 2 + C 1 ϵ 0 2 2 ) .
As a consequence of Lemma 7, we obtain that
ω n 2 ω 0 + t n β g Γ ( 1 + β ) E β ( 2 C 1 t n β ) , n I ¯ N ,
for sufficiently small values of τ , which is what we wanted to prove.    □
The uniqueness of numerical solutions readily follows from this theorem under the same hypotheses. More precisely, the next result holds.
Corollary 3.
Let F C 2 ( R ) and F , F L ( R ) . For sufficiently small values of τ satisfying the inequality (51), the numerical model (42) is uniquely solvable, for any set of initial conditions.
The property of convergence of the numerical method can be established in a similar fashion as the property of nonlinear stability, using the bound for the local truncation error stated in Theorem 4. More precisely, the numerical method (42) is convergent with order of convergence O ( τ 3 β + h 2 2 ) . We omit the proof to avoid redundancy.
We will devote the remainder of this section to produce some simulations using a computer implementation of the finite-difference method (42). As we pointed out previously, the numerical model is nonlinear, which means that its implementation will be carried out using the Newton–Raphson technique. In each iteration, we will use a tolerance equal to 1 × 10 6 under the Euclidean norm, and a maximum number of iterations equal to 20. Here, we would like to mention that this number of iterations was never exceeded in our simulations and, in most of the cases, the number of iterations required was less than 10. For the sake of simplicity, we set p = 1 for the remainder of this work. For the sake of convenience, Figure 1 provides the list of Matlab commands used to produce the three-dimensional figures of in this work. It is obvious that the implementation relies on the use of the command surf, which is used to produce three-dimensional surface plots.
The first two examples are qualitative computational experiments. In the first of them, we compare the results of our numerical simulations against known exact analytical solutions of the system (13) in the non-fractional case.
Example 1.
In our first example, we will check the validity of our computer implementation against a known exact solution of the mathematical model (13). To that end, let us consider the case when α = 2 , β = 1.999 and γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . In this case, the resulting model is approximately the classical sine-Gordon equation without damping, a model for which the following function is an exact static breather solution in the ( 1 + 1 ) -dimensional scenario [58] when the frequency ω R satisfied | ω | < 1 :
ϕ ( x , t ) = 4 arctan 1 ω 2 sin ( ω t ) w cosh ( 1 ω 2 x ) , ( x , t ) R × [ 0 , ) .
We will use this solution to prescribe the initial conditions for the continuous model (13). Computationally, let Ω = ( 15 , 15 ) R , and fix a temporal period T = 20 . The computational parameters are h = 0.5 and τ = 0.005 . Under these circumstances, Figure 2 depicts (a) the numerical approximation obtained using the discrete model (42) and (b) the exact solution for the classical sine-Gordon equation without damping. It is worth noticing that the numerical results are in very good agreement with the exact analytical solution. For convenience, Figure 3 presents graphs of the absolute error between the exact analytical solution and the numerical results. Notice that the results show that the error in the infinity norm is less than 6 × 10 4 . We conclude that the simulations are in good agreement with respect to the analytical solution, considering that the computational parameter h is relatively large.
The purpose of the next computational experiment is to show qualitatively that the discrete model is capable of preserving the conserved quantity.
Example 2.
In a second example, we will let γ = 0 and F ( u ) = 1 cos ( u ) , and consider various possible values of the parameters α and β. As a starting condition, we employ again the exact solution (92) to prescribe exactly the initial profile and velocity. Figure 4, Figure 5 and Figure 6 thus provide graphs of the numerical approximation to the exact form of Φ versus x and t (left columns), together with plots of the associated functional E versus t (right column), for various pairs of parameter values α and β in ( 1 , 2 ) . Notice that the values associated with the approximations to E are approximately equal to zero, in agreement with the theoretical results derived in this work.
Finally, our last experiment is a set of quantitative comparisons aimed at showing that the finite-difference scheme (42) has order of convergence equal to O ( τ 3 β + h 2 2 ) .
Example 3.
In this last experiment, we will consider the system (13) with Ω = ( 15 , 15 ) , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Let us fix α = 1.5 , and use the function (92) to prescribe exactly the conditions at the initial time. As is usual, we will impose homogeneous Neumann boundary data. Various values of β ( 1 , 2 ) will be considered. Beforehand, we must point out that there is no exact solution to compare against in these cases. For that reason, as a close approximation to the exact solution, we employ the numerical solution obtained through our method for h = 1 × 10 3 and τ = 2 × 10 6 . We focus our attention on the numerical study of temporal convergence in this example. To that end, let ϕ N and Φ N represent, respectively, the exact and the numerical solutions at the time T. We define the absolute error at the time T as
ε τ = ϕ N Φ N 2 .
We also evaluate the standard error rate
ρ τ = log 2 ε 2 τ ε τ .
With these conventions, Table 1 provides the numerical study of convergence for two different values of β ( 1 , 2 ) , and different values for the computational parameters τ and h. The results confirm the fact that the finite-difference scheme (42) has temporal convergence rate approximately equal to τ 3 β . This is in obvious agreement with the theoretical results. To avoid redundancy, the numerical study of spatial convergence is not presented here, but the results showed that the scheme is convergent, with spatial order of convergence equal to h 2 .

5. Conclusions

In this work, we considered a generalized nonlinear wave equation in multiple dimensions which extends well-known conservative models of mathematical physics from their classical to their fractional representation. More precisely, the system considered here is a hyperbolic model which includes fractional temporal derivatives of the Caputo type and spatial Riesz-fractional derivatives. Using the energy method, we first showed that our fractional model possesses the expected energy conservation law. Motivated by this fact, we proposed a finite-difference discretization of the equations of motion in which the Caputo temporal derivative is approximated using the L 1 -scheme, and the Riesz spatial derivative is estimated using fractional-order centered differences. We proved the existence and uniqueness of solutions for our model using a fixed-point theorem, and verified that a consistent discretization of the conserved quantity is also preserved throughout time.
From the numerical point of view, we established theoretically the consistency, the stability, and convergence of the proposed scheme. To establish the last two properties, we used a discrete version of Gronwall’s inequality and some a priori error estimates. The numerical model was implemented computationally, and various simulations were obtained. As a result of our simulations, we established the fact that the proposed discrete quantity associated to the numerical model is conserved throughout time, as predicted by the theoretical findings. Finally, it is worth pointing out that the present manuscript is one of the first works in which a conserved quantity is identified for a Caputo–Riesz time-space-fractional nonlinear hyperbolic system, and in which a discretization with similar properties is provided and analyzed using finite differences.
As one of the reviewers pointed out, there are some other works about numerical methods for fractional wave models. Their methods are based on finite differences [54], finite elements [59], spectral techniques [60], among other approaches [61,62,63,64]. It is worth mentioning, however, that these studies present some limitations with respect to the present work. In fact, some of them only consider fractional derivatives in space or time, but not both, while others establish the stability and convergence of the schemes but neglect to study the variational properties of the discretizations. On the other hand, some papers in the literature study the energy conservation of the schemes, but neglect the investigation of the numerical properties of their methods. In fact, the present work is a complete study of space-time-fractional hyperbolic model in which both energy and numerical properties are rigorously established.
One of the most important and novel results in this work is that we prove mathematically the stability and convergence for our finite-difference schemes, while previous reports in the literature only investigate parabolic models. However, in the present manuscript, the investigation was carried out for the hyperbolic case. This was perhaps the most important difficulty encountered by the authors. The way we bypassed this obstacle was by employing a suitable discrete form of Gronwall’s inequality, previously used in the parabolic case. In the present paper, we were able to apply it conveniently in the hyperbolic scenario, in order to prove the stability and the convergence of our discrete model. Moreover, all numerical simulations presented in this work confirmed the validity of the theoretical results, including the conservation law and the convergence rate.

Author Contributions

Conceptualization, J.E.M.-D.; data curation, J.E.M.-D.; formal analysis, J.E.M.-D.; funding acquisition, J.E.M.-D. and T.B.; investigation, J.E.M.-D.; methodology, J.E.M.-D.; project administration, J.E.M.-D.; resources, J.E.M.-D.; software, J.E.M.-D.; supervision, J.E.M.-D. and T.B.; validation, J.E.M.-D. and T.B.; visualization, J.E.M.-D.; writing—original draft, J.E.M.-D.; writing—review and editing, J.E.M.-D., T.B. All authors have read and agreed to the published version of the manuscript.

Funding

J.E.M.-D.: 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 for Science and Technology of Mexico (CONACYT) through grant A1-S-45928. T.B.: Financial support for this work was provided under the scientific project No. 21-71-30011 of the Russian Science Foundation.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors wish to thank the anonymous reviewers for their comments and questions. All of their suggestions resulted in a substantial improvement of the quality of our paper. One of us (J.E.M.-D.) wishes to thank A. J. Serna-Reyes for pointing out various typos and errors in a preliminary version of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Singh, J.; Kumar, D.; Hammouch, Z.; Atangana, A. A fractional epidemiological model for computer viruses pertaining to a new fractional derivative. Appl. Math. Comput. 2018, 316, 504–515. [Google Scholar] [CrossRef]
  2. Qureshi, S. Real life application of Caputo fractional derivative for measles epidemiological autonomous dynamical system. Chaos Solitons Fractals 2020, 134, 109744. [Google Scholar] [CrossRef]
  3. Qureshi, S.; Yusuf, A.; Shaikh, A.A.; Inc, M. Transmission dynamics of varicella zoster virus modeled by classical and novel fractional operators using real statistical data. Phys. A Stat. Mech. Its Appl. 2019, 534, 122149. [Google Scholar] [CrossRef]
  4. Barros, L.C.d.; Lopes, M.M.; Pedro, F.S.; Esmi, E.; Santos, J.P.C.d.; Sánchez, D.E. The memory effect on fractional calculus: An application in the spread of COVID-19. Comput. Appl. Math. 2021, 40, 1–21. [Google Scholar] [CrossRef]
  5. Alshomrani, A.S.; Ullah, M.Z.; Baleanu, D. Caputo SIR model for COVID-19 under optimized fractional order. Adv. Differ. Equ. 2021, 2021, 185. [Google Scholar] [CrossRef] [PubMed]
  6. Ghanbari, B.; Atangana, A. A new application of fractional Atangana–Baleanu derivatives: Designing ABC-fractional masks in image processing. Phys. A Stat. Mech. Its Appl. 2020, 542, 123516. [Google Scholar] [CrossRef]
  7. Qureshi, S.; Yusuf, A.; Shaikh, A.A.; Inc, M.; Baleanu, D. Fractional modeling of blood ethanol concentration system with real data application. Chaos Interdiscip. J. Nonlinear Sci. 2019, 29, 013143. [Google Scholar] [CrossRef]
  8. Ghanbari, B.; Günerhan, H.; Srivastava, H. An application of the Atangana-Baleanu fractional derivative in mathematical biology: A three-species predator-prey model. Chaos Solitons Fractals 2020, 138, 109910. [Google Scholar] [CrossRef]
  9. Ming, H.; Wang, J.; Fečkan, M. The application of fractional calculus in Chinese economic growth models. Mathematics 2019, 7, 665. [Google Scholar] [CrossRef]
  10. Hassani, H.; Machado, J.T.; Mehrabi, S. An optimization technique for solving a class of nonlinear fractional optimal control problems: Application in cancer treatment. Appl. Math. Model. 2021, 93, 868–884. [Google Scholar] [CrossRef]
  11. Ortigueira, M.D. A New Look at the Initial Condition Problem. Mathematics 2022, 10, 1771. [Google Scholar] [CrossRef]
  12. Ortigueira, M.D. An Entropy Paradox Free Fractional Diffusion Equation. Fractal Fract. 2021, 5, 236. [Google Scholar] [CrossRef]
  13. Yavuz, M.; Özdemir, N. Comparing the new fractional derivative operators involving exponential and Mittag-Leffler kernel. Discret. Contin. Dyn. Syst.-S 2020, 13, 995. [Google Scholar] [CrossRef]
  14. Saad, K.M.; Baleanu, D.; Atangana, A. New fractional derivatives applied to the Korteweg–de Vries and Korteweg–de Vries–Burger’s equations. Comput. Appl. Math. 2018, 37, 5203–5216. [Google Scholar] [CrossRef]
  15. Ortigueira, M.D. Two-sided and regularised Riesz-Feller derivatives. Math. Methods Appl. Sci. 2021, 44, 8057–8069. [Google Scholar] [CrossRef]
  16. Ortigueira, M.D.; Bengochea, G. Bilateral tempered fractional derivatives. Symmetry 2021, 13, 823. [Google Scholar] [CrossRef]
  17. Muslih, S.I.; Agrawal, O.P. Riesz fractional derivatives and fractional dimensional space. Int. J. Theor. Phys. 2010, 49, 270–275. [Google Scholar] [CrossRef]
  18. Tarasov, V.E.; Zaslavsky, G.M. Fractional dynamics of systems with long-range interaction. Commun. Nonlinear Sci. Numer. Simul. 2006, 11, 885–898. [Google Scholar] [CrossRef]
  19. Tarasov, V.E. Continuous limit of discrete systems with long-range interaction. J. Phys. A Math. Gen. 2006, 39, 14895. [Google Scholar] [CrossRef]
  20. Christodoulidi, H.; Bountis, A.; Drossos, L. The effect of long-range interactions on the dynamics and statistics of 1D Hamiltonian lattices with on-site potential. Eur. Phys. J. Spec. Top. 2018, 227, 563–573. [Google Scholar] [CrossRef]
  21. Macías-Díaz, J.E.; Bountis, A. Supratransmission in β-Fermi–Pasta–Ulam chains with different ranges of interactions. Commun. Nonlinear Sci. Numer. Simul. 2018, 63, 307–321. [Google Scholar] [CrossRef]
  22. Macías-Díaz, J.E.; Bountis, A. Nonlinear supratransmission in quartic Hamiltonian lattices with globally interacting particles and on-site potentials. J. Comput. Nonlinear Dyn. 2021, 16, 021001. [Google Scholar] [CrossRef]
  23. Caputo, M. The role of memory in modeling social and economic cycles of extreme events. In A Handbook of Alternative Theories of Public Economics; Edward Elgar Publishing: Cheltenham, UK, 2014; pp. 245–259. [Google Scholar]
  24. Tarasova, V.V.; Tarasov, V.E. Elasticity for economic processes with memory: Fractional differential calculus approach. Fract. Differ. Calc. 2016, 6, 219–232. [Google Scholar] [CrossRef]
  25. Škovránek, T.; Podlubny, I.; Petráš, I. Modeling of the national economies in state-space: A fractional calculus approach. Econ. Model. 2012, 29, 1322–1327. [Google Scholar] [CrossRef]
  26. Jiang, H.; Liu, F.; Turner, I.; Burrage, K. Analytical solutions for the multi-term time–space Caputo–Riesz fractional advection–diffusion equations on a finite domain. J. Math. Anal. Appl. 2012, 389, 1117–1127. [Google Scholar] [CrossRef]
  27. Chen, M.; Deng, W.; Wu, Y. Superlinearly convergent algorithms for the two-dimensional space–time Caputo–Riesz fractional diffusion equation. Appl. Numer. Math. 2013, 70, 22–41. [Google Scholar] [CrossRef]
  28. Shen, S.; Liu, F.; Anh, V. Numerical approximations and solution techniques for the space-time Riesz–Caputo fractional advection-diffusion equation. Numer. Algorithms 2011, 56, 383–403. [Google Scholar] [CrossRef]
  29. Shen, Y.; El-Dib, Y.O. A periodic solution of the fractional sine-Gordon equation arising in architectural engineering. J. Low Freq. Noise Vib. Act. Control 2021, 40, 683–691. [Google Scholar] [CrossRef]
  30. Bernard, D.; Leclair, A. The fractional supersymmetric sine-Gordon models. Phys. Lett. B 1990, 247, 309–316. [Google Scholar] [CrossRef]
  31. Altybay, A.; Ruzhansky, M.; Sebih, M.E.; Tokmagambetov, N. Fractional Klein-Gordon equation with singular mass. Chaos Solitons Fractals 2021, 143, 110579. [Google Scholar] [CrossRef]
  32. Laskin, N. Fractional schrödinger equation. Phys. Rev. E 2002, 66, 056108. [Google Scholar] [CrossRef] [PubMed]
  33. Macías-Díaz, J.E. Existence of solutions of an explicit energy-conserving scheme for a fractional Klein–Gordon–Zakharov system. Appl. Numer. Math. 2020, 151, 40–43. [Google Scholar] [CrossRef]
  34. Jin, B.; Lazarov, R.; Zhou, Z. An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. 2016, 36, 197–221. [Google Scholar] [CrossRef]
  35. Ortigueira, M.D. Fractional central differences and derivatives. IFAC Proc. Vol. 2006, 39, 58–63. [Google Scholar] [CrossRef]
  36. Li, M.; Huang, C.; Zhao, Y. Fast conservative numerical algorithm for the coupled fractional Klein-Gordon-Schrödinger equation. Numer. Algorithms 2020, 84, 1081–1119. [Google Scholar] [CrossRef]
  37. Wang, J.; Xiao, A. Conservative Fourier spectral method and numerical investigation of space fractional Klein–Gordon–Schrödinger equations. Appl. Math. Comput. 2019, 350, 348–365. [Google Scholar] [CrossRef]
  38. Wang, P.; Huang, C. An energy conservative difference scheme for the nonlinear fractional Schrödinger equations. J. Comput. Phys. 2015, 293, 238–251. [Google Scholar] [CrossRef]
  39. Duo, S.; Zhang, Y. Mass-conservative Fourier spectral methods for solving the fractional nonlinear Schrödinger equation. Comput. Math. Appl. 2016, 71, 2257–2271. [Google Scholar] [CrossRef]
  40. Martínez, R.; Macías-Díaz, J.E. An energy-preserving and efficient scheme for a double-fractional conservative Klein–Gordon–Zakharov system. Appl. Numer. Math. 2020, 158, 292–313. [Google Scholar] [CrossRef]
  41. Serna-Reyes, A.J.; Macías-Díaz, J.E. Theoretical analysis of a conservative finite-difference scheme to solve a Riesz space-fractional Gross–Pitaevskii system. J. Comput. Appl. Math. 2022, 404, 113413. [Google Scholar] [CrossRef]
  42. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998; Volume 198. [Google Scholar]
  43. Alikhanov, A. A priori estimates for solutions of boundary value problems for fractional-order equations. Differ. Equ. 2010, 46, 660–666. [Google Scholar] [CrossRef]
  44. Barone, A.; Esposito, F.; Magee, C.; Scott, A. Theory and applications of the sine-Gordon equation. La Riv. Del Nuovo Cimento (1971–1977) 1971, 1, 227–267. [Google Scholar] [CrossRef]
  45. Nunes, R.S.O.; Bastos, W.D. Energy decay for the linear Klein–Gordon equation and boundary control. J. Math. Anal. Appl. 2014, 414, 934–944. [Google Scholar] [CrossRef]
  46. Ginibre, J.; Velo, G. The global Cauchy problem for the non linear Klein–Gordon equation. Math. Z. 1985, 189, 487–505. [Google Scholar] [CrossRef]
  47. Popov, C.A. Perturbation theory for the double sine-Gordon equation. Wave Motion 2005, 42, 309–316. [Google Scholar] [CrossRef]
  48. Macías-Díaz, J.E.; Medina-Ramírez, I.E. An implicit four-step computational method in the study on the effects of damping in a modified α-Fermi–Pasta–Ulam medium. Commun. Nonlinear Sci. Numer. Simul. 2009, 14, 3200–3212. [Google Scholar] [CrossRef]
  49. Zhuang, P.; Liu, F. Implicit difference approximation for the time fractional diffusion equation. J. Appl. Math. Comput. 2006, 22, 87–99. [Google Scholar] [CrossRef]
  50. Lin, Y.; Xu, C. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. Phys. 2007, 225, 1533–1552. [Google Scholar] [CrossRef]
  51. Ortigueira, M.D. Riesz potential operators and inverses via fractional centred derivatives. Int. J. Math. Math. Sci. 2006, 2006, 048391. [Google Scholar] [CrossRef]
  52. Çelik, C.; Duman, M. Crank–Nicolson method for the fractional diffusion equation with the Riesz fractional derivative. J. Comput. Phys. 2012, 231, 1743–1750. [Google Scholar] [CrossRef]
  53. Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef]
  54. Macías-Díaz, J.E. On the solution of a Riesz space-fractional nonlinear wave equation through an efficient and energy-invariant scheme. Int. J. Comput. Math. 2019, 96, 337–361. [Google Scholar] [CrossRef]
  55. Browder, F. Existence and uniqueness theorems for solutions of nonlinear boundary value problems. Appl. Nonlinear Partial Differ. Eq. Inn Math. Phys. 1965, 17, 24–49. [Google Scholar]
  56. Sun, Z.z.; Wu, X. A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math. 2006, 56, 193–209. [Google Scholar] [CrossRef]
  57. Hendy, A.S.; Macías-Díaz, J.E. A discrete Grönwall inequality and energy estimates in the analysis of a discrete model for a nonlinear time-fractional heat equation. Mathematics 2020, 8, 1539. [Google Scholar] [CrossRef]
  58. Ablowitz, M.J.; Kaup, D.J.; Newell, A.C.; Segur, H. Method for solving the sine-Gordon equation. Phys. Rev. Lett. 1973, 30, 1262. [Google Scholar] [CrossRef]
  59. Wang, J.; Yin, B.; Liu, Y.; Li, H.; Fang, Z. Mixed finite element algorithm for a nonlinear time fractional wave model. Math. Comput. Simul. 2021, 188, 60–76. [Google Scholar] [CrossRef]
  60. Wang, N.; Shi, D. Two efficient spectral methods for the nonlinear fractional wave equation in unbounded domain. Math. Comput. Simul. 2021, 185, 696–718. [Google Scholar] [CrossRef]
  61. Pandit, S.; Mittal, R. A numerical algorithm based on scale-3 Haar wavelets for fractional advection dispersion equation. Eng. Comput. 2020, 38, 1706–1724. [Google Scholar] [CrossRef]
  62. Mittal, R.; Pandit, S. A Numerical Algorithm to Capture Spin Patterns of Fractional Bloch Nuclear Magnetic Resonance Flow Models. J. Comput. Nonlinear Dyn. 2019, 14, 081001. [Google Scholar] [CrossRef]
  63. Mittal, R.; Pandit, S. Quasilinearized Scale-3 Haar wavelets-based algorithm for numerical simulation of fractional dynamical systems. Eng. Comput. 2018, 35, 1907–1931. [Google Scholar] [CrossRef]
  64. Shukla, H.; Tamsir, M.; Jiwari, R.; Srivastava, V.K. A numerical algorithm for computation modelling of 3D nonlinear wave equations based on exponential modified cubic B-spline differential quadrature method. Int. J. Comput. Math. 2018, 95, 752–766. [Google Scholar] [CrossRef]
Figure 1. Function used to obtain the three-dimensional figures in our numerical experiments. The function uses Matlab commands, and the results were obtained using this computer software.
Figure 1. Function used to obtain the three-dimensional figures in our numerical experiments. The function uses Matlab commands, and the results were obtained using this computer software.
Fractalfract 06 00500 g001
Figure 2. Comparison of the numerical solution obtained using the finite-difference method (42) against the exact solution of the mathematical model (13). In the continuous model, we considered the parameters α = β = 2 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . The exact solution is given by the function (92). Meanwhile, for the numerical approximations we used the same model parameters with β = 1.99 , Ω = ( 15 , 15 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 . For convenience, the first row presents the three-dimensional surfaces for the numerical solution (left column) and the exact solution (right column). Meanwhile, the second row shows the respective interpolated checkerboard plots.
Figure 2. Comparison of the numerical solution obtained using the finite-difference method (42) against the exact solution of the mathematical model (13). In the continuous model, we considered the parameters α = β = 2 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . The exact solution is given by the function (92). Meanwhile, for the numerical approximations we used the same model parameters with β = 1.99 , Ω = ( 15 , 15 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 . For convenience, the first row presents the three-dimensional surfaces for the numerical solution (left column) and the exact solution (right column). Meanwhile, the second row shows the respective interpolated checkerboard plots.
Fractalfract 06 00500 g002
Figure 3. Graph of the point-wise absolute error between the exact analytical solution obtained in Figure 2 and the corresponding numerical approximation. The error was calculated in the infinity norm of R 2 at each point on the spatial domain and at the time T = 20 . For convenience, we provide the three-dimensional plot (left column) and the corresponding interpolated checkerboard plot.
Figure 3. Graph of the point-wise absolute error between the exact analytical solution obtained in Figure 2 and the corresponding numerical approximation. The error was calculated in the infinity norm of R 2 at each point on the spatial domain and at the time T = 20 . For convenience, we provide the three-dimensional plot (left column) and the corresponding interpolated checkerboard plot.
Fractalfract 06 00500 g003
Figure 4. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.2 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Figure 4. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.2 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Fractalfract 06 00500 g004
Figure 5. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.6 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Figure 5. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.6 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Fractalfract 06 00500 g005
Figure 6. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.8 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Figure 6. Left column: graphs of the numerical solution for the mathematical model (13). Right column: graphs of the approximate dynamics for the associated functional E . The parameters employed were α = 1.8 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . Moreover, we used β = 1.2 (top row), β = 1.6 (middle row) and β = 1.8 (bottom row). We employed the finite-difference method (42) with Ω = ( 20 , 20 ) and T = 20 . Computationally, we let h = 0.5 and τ = 0.005 .
Fractalfract 06 00500 g006
Table 1. Numerical study of the convergence of the finite-difference method (42). In this experiment, we let Ω = ( 15 , 15 ) , α = 1.5 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . We used homogeneous Dirichlet boundary conditions, and the initial data were prescribed by the function (92) around the time t = 0 . Two values of β were used, namely, (a) β = 1.2 and (b) β = 1.8 . Various values for the computational parameters τ and h were considered. The results confirm that the numerical method (42) has temporal order of convergence equal to τ 3 β , in agreement with the theoretical results derived in this work.
Table 1. Numerical study of the convergence of the finite-difference method (42). In this experiment, we let Ω = ( 15 , 15 ) , α = 1.5 , γ = 0 and F ( u ) = 1 cos ( u ) , for each u R . We used homogeneous Dirichlet boundary conditions, and the initial data were prescribed by the function (92) around the time t = 0 . Two values of β were used, namely, (a) β = 1.2 and (b) β = 1.8 . Various values for the computational parameters τ and h were considered. The results confirm that the numerical method (42) has temporal order of convergence equal to τ 3 β , in agreement with the theoretical results derived in this work.
(a) β = 1.2
h = 0 . 5 h = 0 . 25
τ ϵ τ ρ τ ϵ τ ρ τ
1 × 10 3 2.478298 × 10 4 6.667961 × 10 5
5 × 10 4 7.739979 × 10 5 1.492801 2.216056 × 10 5 1.589251
2.5 × 10 4 2.574133 × 10 5 1.588243 6.391346 × 10 6 1.793803
1.25 × 10 4 8.276150 × 10 6 1.637055 1.781890 × 10 6 1.842711
6.25 × 10 5 2.391394 × 10 6 1.791108 5.408490 × 10 7 1.720111
(b) β = 1 . 8
h = 0 . 5 h = 0 . 25
τ ϵ τ ρ τ ϵ τ ρ τ
5 × 10 4 9.279493 × 10 5 2.202533 × 10 5
2.5 × 10 4 5.031979 × 10 5 0.882920 1.071691 × 10 5 1.039274
1.25 × 10 4 2.581205 × 10 5 0.963081 4.764840 × 10 6 1.169390
6.25 × 10 5 1.259219 × 10 5 1.035515 2.006268 × 10 6 1.247913
3.125 × 10 5 5.742301 × 10 6 1.132829 8.809504 × 10 7 1.187382
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Macías-Díaz, J.E.; Bountis, T. An Efficient Dissipation-Preserving Numerical Scheme to Solve a Caputo–Riesz Time-Space-Fractional Nonlinear Wave Equation. Fractal Fract. 2022, 6, 500. https://doi.org/10.3390/fractalfract6090500

AMA Style

Macías-Díaz JE, Bountis T. An Efficient Dissipation-Preserving Numerical Scheme to Solve a Caputo–Riesz Time-Space-Fractional Nonlinear Wave Equation. Fractal and Fractional. 2022; 6(9):500. https://doi.org/10.3390/fractalfract6090500

Chicago/Turabian Style

Macías-Díaz, Jorge E., and Tassos Bountis. 2022. "An Efficient Dissipation-Preserving Numerical Scheme to Solve a Caputo–Riesz Time-Space-Fractional Nonlinear Wave Equation" Fractal and Fractional 6, no. 9: 500. https://doi.org/10.3390/fractalfract6090500

APA Style

Macías-Díaz, J. E., & Bountis, T. (2022). An Efficient Dissipation-Preserving Numerical Scheme to Solve a Caputo–Riesz Time-Space-Fractional Nonlinear Wave Equation. Fractal and Fractional, 6(9), 500. https://doi.org/10.3390/fractalfract6090500

Article Metrics

Back to TopTop