Numerically Efﬁcient Methods for Variational Fractional Wave Equations: An Explicit Four-Step Scheme

: In this work, we investigate numerically a one-dimensional wave equation in generalized form. The system considers the presence of constant damping and functional anomalous diffusion of the Riesz type. Reaction terms are also considered, in such way that the mathematical model can be presented in variational form when damping is not present. As opposed to previous efforts available in the literature, the reaction terms are not only functions of the solution. Instead, we consider the presence of smooth functions that depend on fractional derivatives of the solution function. Using a ﬁnite-difference approach, we propose a numerical scheme to approximate the solutions of the fractional wave equation. Along with this integrator, we propose discrete forms of the local and the total energy operators. In a ﬁrst stage, we show rigorously that the energy properties of the continuous system are mimicked by our discrete methodology. In particular, we prove that the discrete system is dissipative (respectively, conservative) when damping is present (respectively, absent), in agreement with the continuous model. The theoretical numerical analysis of this system is more complicated in light of the presence of the functional form of the anomalous diffusion. To solve this problem, some novel technical lemmas are proved and used to establish the stability and the quadratic convergence of the scheme. Finally, we provide some computer simulations to show the capability of the scheme to conserve/dissipate the energy. Various fractional problems with functional forms of the anomalous diffusion of the solution are considered to that effect.


Introduction
The design of energy-preserving methods for physical systems has been a fruitful avenue of research in the last decades. Historically, the problem of designing energy-conserving methods may date back to the decade of the 1970s [1,2]. However, L. Vázquez and coauthors were probably the first researchers who pointed out the physical and mathematical significance of designing this type of schemes [3]. Various seminal papers by Vázquez and his coworkers were published in the 1990s, including various energy-conserving numerical schemes to solve partial differential equations such as the Schrödinger equation [4], the sine-Gordon equation [5,6], the Klein-Gordon equations [7], and even systems consisting of ordinary differential equations [8]. In those papers, the authors established thoroughly the capability of their schemes to preserve the energy properties of the continuous problem. Moreover, they employed a discrete form of the energy method to establish rigorously the stability and the convergence properties of the schemes. After the publication of those works, the investigation and speed, we implemented the numerical method in a Fortran 95. However, we must point out that the algorithm can be implemented in any scientific computer language, including C++ or Matlab.
This work is organized as follows. In Section 2, we present the fractional partial differential equation of interest. The model is an extension of some continuous generalizations of the classical wave equation with constant damping and fractional diffusion. We propose an energy density function associated to our model. Moreover, we prove therein that the total energy of the system is conserved or dissipated, depending on whether damping is absent or present. Section 3 introduces the discrete nomenclature employed throughout this work, and presents the finite-difference scheme to approximate the solutions of the mathematical model. Discrete forms of the local and total energy are presented therein, and we establish a discrete analog of the theorem on the conservation/dissipation of energy of the continuous model. As expected, the boundedness of the numerical solutions is established under the same conditions of the continuous-case scenario. In turn, the purpose of Section 4 is to establish the main numerical features of our discrete model. Some illustrative simulations are shown in Section 5. We close this work with some concluding remarks.

Preliminaries
We let I n = {1, . . . , n} and I n = I n ∪ {0}, for each n ∈ N. Throughout, we let γ ∈ R + ∪ {0} represent physically a constant damping coefficient. Let a and b be real numbers such that a < b, and let B = (a, b) ⊆ R and Ω = B × (0, T) ⊆ R 2 . We employ the symbols ∂B, B and Ω to denote, respectively, the boundary of B, and the closures of B and Ω with respect to the standard topology of R 2 . In the following, we let u : Ω → R be a function, and extend the definition of u by letting u(x, t) = 0, for each (x, t) ∈ (R \ B) × [0, T]. Throughout this manuscript, we let Γ denote the usual Gamma function that generalizes factorials.

Definition 1.
Let L x,2 (Ω) denote the set of all functions f : Ω → R such that f (·, t) ∈ L 2 (B) for each t ∈ [0, T]. Moreover, for each pair f , g ∈ L x,2 (Ω), the inner product of f and g is the function of t defined by In turn, the Euclidean norm of f ∈ L x,2 (Ω) is the function of t defined by f x,2 = f , f x . In general, if 1 ≤ q < ∞, then L x,q (Ω) represents the set of all functions f : Ω → R such that f (·, t) ∈ L q (B), for each t ∈ [0, T]. For each such function f , we define its norm as the function of t given by Definition 2. Let f : R → R be a function, and let n ∈ N ∪ {0} and α ∈ R satisfy n − 1 < α < n. The Riesz fractional derivative of f of order α at x ∈ R is defined (when it exists) as Definition 3. Let α ∈ R and n ∈ N be such that n − 1 < α < n. When it exists, the Riesz space-fractional partial derivative of the function u : Ω → R of order α with respect to x at the point (x, t) is defined by For the sake of convenience, if n ∈ {1, 2}, then we let the Riesz partial derivative of u of order n to be equal to the usual partial derivative of u of order n with respect to x. Riesz fractional partial derivatives are self-adjoint and non-positive operators [48]. It follows that their additive inverses are self-adjoint and nonnegative, which implies, in turn, that they possess square-root operators [49]. This means in particular that, if α ∈ (1, 2] and u, v : Ω → R, then For the remainder of this manuscript, we let φ, ψ : Ω → R be sufficiently smooth functions vanishing at the boundary of Ω, and let p ∈ N. Suppose that α 1 , . . . , α p ∈ (1, 2] are constants, and let G, G 1 , G 2 , . . . , G p : R → R be smooth functions. In particular, we assume that G = G(v) and G j = G j (v), for each v ∈ R and j ∈ I p . With these conventions, the problem under investigation in this work is the initial-value problem with homogeneous Dirichlet boundary data Definition 4. Let u be a function satisfying the initial-boundary-value problem in Equation (6). We define the local energy density of the undamped system as the function In turn, the total energy at the time t ∈ (0, T) is given by Theorem 1 (Dissipation of energy). If u satisfies Equation (6), then E (t) = −γ u t 2 x,2 , where u t = u t (x, t) is the partial derivative of u with respect to t at (x, t). The system in Equation (6) is conservative if γ = 0. In general, Proof. We take the derivative of each of the three terms at the right-hand side of Equation (8). For the third term, we use the square-root properties of fractional derivatives along with the chain rule. The first two terms only require the use of the chain rule. In such way, we readily obtain and for each t ∈ (0, T) and j ∈ I p . Taking the derivative with respect to t on both sides of Equation (8), exchanging the derivative and integral operators, using then the identities above and regrouping, we reach as desired. It is obvious now that the system is conservative if γ = 0. Moreover, the identity in Equation (9) is also an immediate consequence of the fact that E (t) = −γ u t 2 x,2 , which is obtained by integrating both sides of this equation over (0, t).
Corollary 1 (Energy positivity). Suppose that G, G 1 , G 2 , . . . , G p are all nonnegative functions. Then, E is likewise nonnegative in (0, T) and, moreover, Corollary 2 (Boundedness). If G, G 1 , G 2 , . . . , G p are nonnegative functions, then there exists a constant C ≥ 0 that depends only on φ and ψ, such that each of the terms at the right-hand side of Equation (14) is uniformly bounded by C, for all t ≥ 0.
Proof. Using the hypotheses of this result and Theorem 1, we readily obtain the next inequalities, valid for all t ≥ 0: Let C be the right-hand side of his chain of inequalities. Note that whence the conclusion of this proposition readily follows.
In the following examples, we provide some fractional generalizations of well-known variational systems from the physical sciences. The systems considered therein will be particular forms of the hyperbolic partial differential equation in Equation (6).

Example 1.
In the following examples, we assume that γ = 0 and p = 1.
(a) Suppose that G(v) = 1 − cos(v) and G 1 (v) = 1 2 v 2 , for each v ∈ R. Then, the resulting partial differential equation in Equation (6) is the fractional sine-Gordon equation, which is a well-known physical model that appears in relativistic quantum mechanics. (b) The fractional form of the nonlinear Klein-Gordon equation is obtained from the partial differential equation of Equation (6) when G(v) = 1 2! v 2 − 1 4! v 4 for each v ∈ R, and G 1 is as in (a). The Klein-Gordon equation is also a useful model in relativistic quantum mechanics and particle physics. (c) If G 1 is as in (a) and G(v) = 1 − 1 3 cos v − 1 6 cos(2v) for each v ∈ R, then the partial differential equation resulting in Equation (6) is the fractional double sine-Gordon equation. (d) Let > 0 and p ∈ N satisfy p ≥ 3. If G ≡ 0 and G 1 (v) = 1 2 v 2 + v p for each v ∈ R, then the resulting equation is a fractional and continuous extension of the Fermi-Pasta-Ulam-Tsingou chains [50,51].
Obviously, all the models in this example reduce to their well-known integer-order systems when α = 2.

Numerical Method
In this section, we provide the discrete nomenclature used for the remainder of the present work, and introduce the finite-difference method to solve Equation (6) along with a discrete form of the Hamiltonian in Equation (7). Moreover, we show that a discrete form of Theorem 1 is satisfied in the discrete-case scenario. To that end, we consider a regular partition of the interval B, consisting of M ∈ N subintervals. The norm of the partition is represented by h = (b − a)/M, and the nodes are defined by x m = a + mh, for each m ∈ I M . Similarly, fix a uniform partition of [0, T] consisting of N ∈ N subintervals, and let τ be the respective partition norm. For each n ∈ I N , agree that t n = nτ.
Throughout, we fix J = I M−1 × I N−2 and J = I M × I N . Let us use R h to represent the spatial grid {x m : m ∈ I M }, and V h to represent the real vector space of all functions f : For simplicity, if f ∈ V h and m ∈ I M , then we set f m = f (x m ). In this work, unless we mention something different, we use the symbol U n m to represent an approximation to the exact value of the solution u of Equation (6) at the point (x m , t n ), for each (m, n) ∈ J. Moreover, let The Euclidean norm induced by ·, · is denoted by · 2 , and · ∞ : V h → R is the usual infinity norm in V h , which is defined as U ∞ = max{|U m | : m ∈ I M }, for each U ∈ V h . Definition 6. Let (U n ) n∈I N ⊆ V h . We introduce the linear difference operators on V h defined by For convenience, we let t n+ 1 2 = (n + 1 2 )τ, for each n ∈ I N−1 . If U : Ω → R is a sufficiently smooth function then the value Equation (19) is a second-order consistent approximation for the partial derivative of U with respect to t at the point (x m , t n+ 1 2 ). In addition, the quantity in Equation (20) is a second-order approximation to the value of U at the same point. In turn, Equation (21) approximates the second-order partial derivative of U with respect to t at (x m , t n ) with second order of consistency. It is easy to see thatδ approximates consistently the second-order derivative of U with respect to t at the point (x m , t n+ 1 2 ) with a quadratic order of consistency.
(a) The coefficients (g k . When α ∈ (0, 1) ∪ (1, 2], h > 0 and f ∈ C 5 (R) has its derivatives up to order five which belong to L 1 (R), then the following consistency property holds true (see [52]): Definition 8. Let U ∈ V h and α ∈ (0, 2]. If m = 0, M, then we convey that δ (α) Definition 9. Let (U n ) n∈I N ⊆ V h , and suppose that G : R p+2 → R is a continuously differentiable function. We introduce the nonlinear difference operators δ U G and δ U G j defined on V h by the formulas and for each (m, n) ∈ J and j ∈ I p . In addition, for any function F : R → R, we employ F(U n ) to represent the vector (F(U n 0 ), F(U n 1 ), . . . , F(U n M )) ∈ V h . Moreover, we use the symbol "1" to represent both the multiplicative identity of R and the (M + 1)-dimensional vector all of whose components are equal to 1.
At this stage of our discussion, the nomenclature introduced thus far suffices to provide the full finite-difference discretization of the model in Equation (6). The numerical method proposed in this manuscript is given by the discrete system Here, the functions φ U , ψ U , χ U : B → R are used to prescribe exactly the approximations at the times t 0 , t 1 and t 2 , respectively. The forward-difference stencil of the numerical model in Equation (30) is shown in Figure 1. Obviously, the finite-difference scheme is an explicit four-step technique, whence the existence and the uniqueness of numerical solutions is guaranteed for any set of initial data.
Forward-difference stencil for the approximation to the exact solution of the partial differential equation of Equation (6) at the time t n , using the finite-difference scheme in Equation (30). The black circles represent the known approximations at the times t n−1 , t n and t n+1 , while the cross denotes the unknown approximation at the time t n+2 .
It is important to mention that the approximations at the times t 1 and t 2 can be computed in alternative forms. For example, one may try to approximate the first-order partial derivative with respect to time at the time t = 0 using suitable finite-difference approximations. However, such approach could result in a loss of the quadratic consistency of the numerical model. Of course, this is a shortcoming of the present approach. Nevertheless, the scheme presents many other advantages, which are established thoroughly in the following, including its ability to preserve the energy in the undamped case, its quadratic consistency, its stability and its quadratic convergence. Needless to mention that the numerical model in Equation (30) can be easily implemented in any computer language.
Theorem 2 (Existence and uniqueness of solutions). If G, G 1 , G 2 , . . . , G p ∈ C 1 (R), then the discrete problem in Equation (30) has a unique solution for each set of initial conditions. Definition 10. Let (U n ) n∈I N ⊆ V h be solution of the scheme in Equation (30). We define the discrete energy density as In turn, the total energy at the time t n is given by then the following identities are satisfied for each n ∈ I N−2 : Proof. To establish Identity (a), observe that Identities (b) and (c) are also straightforward, thus we only prove Identity (d). Notice that Identity (d) readily follows now from these identities.
2 , for each n ∈ I N−1 . As a consequence, the discrete system in Equation (30) is conservative if γ = 0. In general, the following identities holds: Proof. We use the fact that (U n ) n∈I N satisfies the scheme in Equation (30) along with the identities of Lemma 3. Collecting terms and substituting the expression of the discrete total energy, we obtain It readily follows that the discrete system is conservative if γ = 0. The first identity of Equation (35) is reached now using Identity (b) of Lemma 3, while the second is obtained from Equation (36).

Numerical Properties
The main numerical properties of the finite-difference scheme in Equation (30) are proved in the present section. More precisely, we establish rigorously the quadratic consistency of the scheme, the stability and its quadratic rate of convergence in both space and time. Various useful results are recalled or proved in the way, including a useful discrete Gronwall inequality and a technical proposition taken from the literature. Moreover, a novel auxiliary lemma is mathematically established, and we employ it to prove the stability and convergence of Equation (30).
In a first stage, we establish the consistency property of our scheme. To that end, we suppose that u : Ω → R is a function, and convey that u n m = u(x m , t n ), for each (m, n) ∈ J. Let us introduce the differential and the difference operators D and D, which are defined by Du n m =δ (2) t u for each (x, t) ∈ Ω and (m, n) ∈ J. In our following results, we observe the discrete nomenclature u n = (u n 0 , u n 1 , . . . , u n M ) ∈ V h , for each n ∈ I M .
Definition 11. If U = (U n ) n∈I N is a sequence in V h , then we define the real number Theorem 4 (Consistency). Let u ∈ C 5,4 x,t (Ω) and suppose that G ∈ C 2 (R) and G 1 , G 2 , . . . , G p ∈ C 5 (R). Then, there exist constants C, C ∈ R + , which are independent of h and τ, such that Proof. Using Taylor's theorem, it is easy to check that there exist constants C 1 , C 2 , C 3 , C 4 ∈ R + , which are independent of h and τ, such that the following hold for each j ∈ I p : 4 : j ∈ I p }. The first inequality of the conclusion follows then from the triangle inequality, letting C = max{C 1 , γC 2 , C 3 , C 4 }. The second inequality is proved analogously.
We turn our attention now to the properties of stability and convergence of the finite-difference scheme in Equation (30). To that end, the following well-known inequalities are needed in the sequel. We use them in the sequel without mentioning them explicitly: (c) If (U n ) N n=0 ⊆ V h and n ∈ I N , then The following auxiliary lemma has been taken from the literature. In the following, we let Lemma 4 (Macías-Díaz [36]). Let G ∈ C 2 (R), and assume that (R n+ 1 2 ) N−1 n=0 is a sequence in V h . Moreover, suppose that (V n ) N n=0 and (W n ) N n=0 are two solutions of Equation (30) corresponding to Φ V and Φ W , respectively. Let ε n = V n − W n for each n ∈ I N , and define Then, there exists a constant C ∈ R + ∪ {0} which depends only on G, such that, for each k ∈ I N−1 , Lemma 5. Suppose that G j ∈ C 2 (R) for each j ∈ I p , and suppose that (V n ) N n=0 and (W n ) N n=0 are solutions of Equation (30) corresponding to Φ V and Φ W , respectively. Let ε n = V n − W n for each n ∈ I N , and define Then, there exists a constant C ∈ R + ∪ {0} that depends only on G, such that, for each k ∈ I N−1 , Proof. Let j ∈ I p . In the following chain of inequalities, we use firstly the square-root properties of fractional centered differences. Then, we employ Lemma 4 to show that there exists a constant C j ∈ R + ∪ {0} that only depends on G j , such that, for each k ∈ I N−1 , Here, C j = C j g (α) h in view of Lemma 2. It is clear that C = C 1 + C 2 + . . . + C p . The second inequality of this result can be established in similar fashion from the second inequality of Lemma 4.
The following discrete version of Gronwall's inequality is of utmost importance.
Lemma 6 (Pen-Yu [54]). Let (ω n ) N n=0 and (ρ n ) N n=0 be finite sequences of nonnegative mesh functions, and suppose that there exists C ≥ 0 such that Then, ω n ≤ ρ n e Cnτ for each n ∈ I N .
Theorem 5 (Stability). Let G, G 1 , G 2 , . . . , G p ∈ C 2 (R). Suppose that Φ V and Φ W are two sets of initial conditions for Equation (30), and that (V n ) n∈I N and (W n ) n∈I N are the respective solutions. Let ε n = V n − W n for each n ∈ I N , and define the nonnegative constants Then, there is a constant C ∈ R + ∪ {0} independent of τ, such that ω k ≤ ρe CT , for each k ∈ I N−1 .
Proof. Notice that the sequences V = (V n ) n∈I N and W = (W n ) n∈I N satisfy the problem in Equation (30) with U = V and U = W, respectively. Subtracting those problems, we obtain the discrete systemδ Let k ∈ I N−2 , and take the inner product of δ t ε n+ 1 2 with the vector difference equation of Equation (57) corresponding to the time t n+ 1 2 . To that end, we apply Identities (a) and (b) of Lemma 3 and rearrange terms. As a consequence, we obtain the expression Take now the sum on both sides of this equation over all indexes n ∈ I k . Use the formula for telescoping sums to simplify algebraically, multiply then by 4τ and add the term 2µ t δ t ε 1 2 2 2 on both sides of the identity. Use then Lemmas 4 and 5 to obtain Here, we let C 1 = C + C , where C and C are the constants in Lemmas 4 and 5, respectively. Now, multiply the nth difference equation of Equation (57) by (−1) n−1 and take the sum for all n ∈ I k . Using then telescoping sums, regrouping and inductions, it is easy to notice that for each ∀k ∈ I N−2 . Solve for δ (2) t ε k+1 and calculate the Euclidean norm of this term. After rearranging terms, simplifying algebraically and using the second inequalities of Lemmas 4 and 5, we obtain for each k ∈ I N−2 . Combining Equations (59) and (61), and using the notation of this theorem, it follows that The conclusion follows now from Lemma 6, letting C = C 1 + C 1 T 3 + 4γ 2 T.
Proof. The proof of this result is similar to that of Theorem 5. For that reason, we only provide an abridged argument. Throughout, let n = u n − U n for each n ∈ I N . Notice that the sequence ( n ) n∈I N satisfies the discrete problem where R ). In light of Theorem 4, it follows that there exists a constant C 0 which is independent of τ, such that | R | ∞ ≤ C 0 (τ 2 + h 2 ). Take the inner product of δ t n+ 1 2 with the nth vector equation of Equation (63), and then sum overall n ∈ I k , for some k ∈ I N−2 . Rearrange terms and use the first inequalities of Lemmas 4 and 5 to obtain that δ t for each k ∈ I N−2 . As in the proof of Theorem 5, we agree that C 1 = C + C , where C and C are the constants in Lemmas 4 and 5, respectively. On the other hand, departing from the difference equation of Equation (63), using mathematical induction and performing some algebraic simplifications lead to the following identity, valid for each ∀k ∈ I N−2 : As a consequence, we notice that Here, C 1 = 5 4 C 1 . Substitute now Equation (66) into Equation (64) and simplify. It is easy to check then that where C 2 = C 1 + C 1 T 3 + 5γ 2 T, and Lemma 6 shows now that ω k ≤ C 3 ρ k for all k ∈ I N−1 , where C 3 = e C 2 T . However, recall that the initial conditions of Equation (63) are zero. In particular, this means that the first three terms on the right-hand side of Equation (69) are equal to zero, for all k ∈ I N−2 . This and the consistency property of the finite-difference scheme in Equation (30) imply that, for each k ∈ I N−1 , Here, C 2 4 = 2C 2 0 C 3 T(4 + 5T). As a consequence, note that Telescoping sums and the fact that 0 = 0, yield the inequality which is what we want to prove. Here, we just need to clarify that C = C 4 T.

Computer Simulations
The purpose of this section is to provide examples that illustrate the validity of Theorem 3. Various different scenarios are considered to that end.

Example 2.
Let us consider the system in Equation (6) with p = 1 and Ω = (0, 125) × (0, 200). We let G ≡ 0 and The resulting model describes a fractional version of a nonlinear vibration problem associated to an elastic string [55]. It is worth pointing out that this problem has interesting applications to the propagation of sounds and mechanical vibrations [56][57][58]. Fix the initial conditions as Computationally, we use the model in Equation (30) with h = 1 and τ = 0.025, and agree that α = α 1 . Under these circumstances, the left column of Figure 2 shows the approximate solutions of the problem for α = 2 (top row), α = 1.6 (middle row) and α = 1.2 (bottom row), using γ = 0. Meanwhile, the right column shows the corresponding discrete energy densities. In turn, Figure 3 shows the dynamics of the total energy of the system for various values of α. These graphs confirm that the energy is conserved in the absence of damping. Figure 4 shows graphs of the total energy of the system for various values of α and γ. The results show that the total energy dissipates when damping is present. These remarks are in agreement with Theorem 3.
It is worth pointing out that our choice of the parameter τ in our previous example was obeying the need to provide a good approximation to the solution. Indeed, the value of τ is much smaller than that of h, but this choice was actually arbitrary. After all, the conditions to guarantee the stability and the convergence of the finite-difference scheme in Equation (30) are independent of the ratio of those parameters.

Example 3.
Let Ω = (0, 128) × (0, 500) and = 0.75, and consider the system in Equation (6) with p = 1, We approximate the exact solutions using Equation (30) with h = 1 and τ = 0.025. As initial conditions, we set φ(x) = v(x, 0), ψ(x) = v(x, τ) and χ(x) = v(x, 2τ), for each x ∈ (0, 128). The function v is the algebraic sum of the well-known kink and anti-kink solutions of the Toda lattice [59], and it is given by the formula v(x, t) = A ln for each (x, t) ∈ R 2 . In our simulations, we let A = 5, κ = 0.1 and α = α 1 . Figure 5 shows the approximate solutions of the corresponding problem in Equation (6) for α = 2 (top row), α = 1.6 (middle row) and α = 1.2 (bottom row), using γ = 0. Note that the dynamics of the solutions and the energy density is similar to that reported in [50]. In turn, Figure 6 shows the dynamics of the total energy of the system, for various values of α and γ. The graphs confirm again that the energy is conserved in the absence of damping, and it is dissipated when damping is present.   The next example confirms the convergence rates derived in Theorem 6. For this purpose, we introduce the maximum-norm error between the exact solution of the continuous problem in Equation (6) at the time T, and the corresponding numerical approximation obtained through Equation (30), which is given by We also define the following standard rates: Example 4. Let us consider the same problem investigated in Example 3, with α = α 1 = 1.5. In particular, we employ the same space-time domain with the same model parameters, the same functions G and G 1 , as well as the same initial data on (0, 128). The exact solution for this problem at the time T = 500 is not known in exact form, but we estimate it using relatively small values of the computational parameters. To that effect, we set h = 4/2 8 and τ = 0.1/2 6 . Under these circumstances, Table 1 shows the maximum-norm errors for various values of τ, keeping the parameter h constant. Additionally, the standard rates in time are provided in each case.
The results show that the method has second order of convergence in time, in agreement with the conclusion of Theorem 6. In turn, Table 2 shows the analysis of spatial convergence of the numerical model. The results confirm the quadratic rate of convergence of Equation (30), in agreement with Theorem 6 again. The following examples provide hard numerical evidence on the quadratic order of convergence of the numerical model in Equation (30), as well as on its stability property.

Example 5.
The purpose of the present example is to provide more solid evidence on the temporal quadratic order of convergence of the scheme in Equation (30). To that end, let p = 1 and Ω = (0, 1) × (0, 1000). Convey that G ≡ 0, and let G 1 be given as in Example 2. As initial data, we set For illustrative purposes, let us set α 1 = 1.6 and γ = 0. Under these conditions, Table 3 provides a temporal analysis of convergence of the finite-difference scheme in Equation (30), considering various values of τ and h. The results of our simulations show that the numerical model exhibits a quadratic order of convergence, even in the case when the parameters h and τ are approximately equal. This result is in full agreement with Theorem 6. Example 6. Recall that one way to establish the stability of variational schemes is to show that the energy of the undamped system remains constant for long periods of time, using different values of the computational parameters. In the present example, we confirm this feature of our numerical model. To that end, we consider the same problem investigated in Example 3. That is, let α = 1.5, γ = 0, B = (0, 128) and = 0.75, and consider the model in Equation (6) with p = 1, G ≡ 0 and G 1 given by Equation (75). 2τ), for each x ∈ (0, 128), where v is provided by Equation (76). Throughout, we set T = 1 × 10 6 , and fix h = 1. Under these circumstances, Figure 7 provides graphs of the dynamics of the total energy of the system versus t, for various values of τ. The range of the scale of the y-axis is the interval [0.17103, 0.17108]. The result confirm that the total energy is approximately constant, as predicted by Theorem 3. These simulations provide support on the fact that the finite-difference scheme in Equation (30) is stable, even for values of τ which are larger than 1, and for all α ∈ (1, 2]. Obviously, this is in agreement with Theorem 5. Before closing this section, it is worth pointing out that we carried out more numerical experiments. We do not include them in this paper to avoid redundancy, but they confirm the validity of the theoretical results obtained in this work. (e) (f) Figure 7. Graphs of the total energy of the system in Equation (6) versus t. The parameters employed are p = 1, G ≡ 0, and G 1 is given by Equation (75). We used = 0.75, α = 1.5, γ = 0, B = (0, 128) and T = 1 × 10 6 . The initial data were defined by the functions φ(x) = v(x, 0), ψ(x) = v(x, τ) and

Conclusions
In this work, we consider a general wave equation that encompasses various mathematical models from the physical sciences, including the sine-Gordon, the Klein-Gordon and the double sine-Gordon equations as well as continuous forms of the classical Fermi-Pasta-Ulam-Tsingou chains. The model under investigation includes the presence of constant damping, a nonlinear reaction term and general functions that depend on anomalous diffusion forms of the solution. In a first step, we show that the undamped form of the system has a variational structure. We propose functions for the local and total energy densities, and we show that the total energy is conserved or dissipated, depending on whether damping is absent or present, respectively. Motivated by these facts, we propose a finite-difference discretization of the continuous model, along with discrete forms of the local and total energy functionals. We show rigorously that the proposed discrete model is capable of reflecting the same energy properties of the continuous model. Moreover, we prove mathematically that the numerical model yields consistent approximations to the exact solutions of the continuous system, with quadratic order of consistency in both space and time.
To establish the properties of stability and convergence, we use various propositions. To start with, a discrete form of Gronwall's inequality available in the literature is employed as the cornerstone of the proofs. Some other technical results are also recalled. However, the proofs of stability and convergence employ a novel technical proposition in order to bound appropriately the functions which depend on anomalous diffusion. As a result, we show that the proposed finite-difference scheme is stable, and that it converges to the exact solution of the continuous model with quadratic order of convergence in both space and time. To prove those results, suitable regularity conditions on the functions involved in the model are needed. For illustration purposes, we provide some computer simulations, which exhibited the capability of the scheme to conserve the energy when damping is absent and dissipate it when damping is present. The computer simulations were obtained using a computer implementation of our finite-difference scheme in Fortran 95. It is worth pointing out that the results of this work could be used in various potential applications, such as the investigation of nonlinear phenomena in generalized fractional systems [53,60], or even the transmission of signals in fractional forms of nonlinear media [61][62][63].
We would like to point out that the present work improves greatly various individual efforts by this author. As mentioned in the Introduction, some methodologies have been proposed already to solve fractional wave equations using fractional-order centered differences [35][36][37].
From the mathematical point of view, those previous efforts considered much simpler forms of the mathematical model in Equation (6). More precisely, in those articles, this author considered the initial-boundary-value problem u(a, t) = u(b, t) = 0, ∀t ∈ (0, T), Notice that the model in the present manuscript is both a dimensional and variational generalization of the problem in Equation (80). Indeed, it is easy to check that the mathematical problem in Equation (80) does not extend the Fermi-Pasta-Ulam-Tsingou problem to the fractional scenario. In that sense, the mathematical model in Equation (6) is a nontrivial extension of those investigated in [35][36][37].
On the other hand, from the numerical point of view, the present methodology presents the advantage of being relatively simple to implement. When compared to the approaches described in our previous papers, the present method is an explicit technique, which has the advantage of being numerically efficient under relatively light conditions. The obvious disadvantage here is the fact that the methodology is a four-step technique, and this fact implies the knowledge of the numerical solution at the first three time-steps. However, on the other hand, this disadvantage is compensated by the fact that the scheme is an energy-preserving method, and that the numerical efficiency analysis results in a relatively straightforward process. Moreover, the conditions under which the scheme is stable and convergent are little demanding in terms of the computational parameters. Additionally, the computer implementation of the scheme is also relatively straightforward, and it can be employed to investigate a wide range of problems in the physical sciences.