Time-space variable-order fractional nonlinear system of thermoelasticity: numerical treatment

This paper focuses on a numerical study of the general time-space variable-order fractional nonlinear problem of thermoelasticity in one dimension using the weighted average nonstandard finite difference (WANSFD). By replacing the second order space derivative with a Riesz fractional variable-order derivative and the time derivative by Caputo fractional variable-order operator in the standard system which arises in thermoelasticity, we obtain this general system. Using a kind of John von Neumann technique, we study the stability of the designed schemes. Also, the truncation error of the introduced schemes is studied. Our numerical treatment is shown graphically. These results expose that WANSFD approach is suitable and effective for solving the proposed system; moreover, it is easy to implement.


Introduction
Today scientists in different fields such as plasma waves, fluid mechanics systems and solid state physics use coupled partial differential equations to describe many phenomena. A system of nonlinear coupled hyperbolic and parabolic equations is always used in studies of circled fuel reactor, radiation hydrodynamics, magnetoelasticity, thermoelasticity, and in biology [1][2][3][4][5]. It is known that studying the behavior of solutions of these systems is a very important and difficult area of research. A lot of authors [4][5][6] have been studying the existence of solutions for such models, as well as their uniqueness and stability.
It is very difficult to obtain the analytic exact solutions of such nonlinear problems and, unfortunately, only in simple cases one can find these solutions. Therefore, it is very important to approximate them numerically, which requires great potential from the researchers. The authors of [7] and of the references therein studied a system of coupled parabolic equations. Some of these studies depended on the finite element technique [8,9] and on the finite difference technique [1, 2, 10] for a nonlinear system of coupled hyperbolic and parabolic equations. In [11,12], the authors apply the Adomian decomposition approach and the variational iteration approach for solving the model. In [3], the authors found the solution of a nonlinear system of coupled hyperbolic and parabolic equations with specified harmonic displacement at the border depending on the Poincaré extension with a small parameter.
It is well known, since the last few years, that the fractional derivative operators are more compatible to describe the properties of many real-life problems and more helpful than those of integer-order to understand such problems [13]. In fact, the next state in the model with a fractional-order derivative is based on all of its prior states and not only on its present state as in the case with a derivative of integer order. And therefore, the derivative of fractional order will be more authentic than that of integer order. Inclusion of the hereditary properties and memory effect of processes in numerous real-life problems is considered as one of the most important advantages of the fractional derivative [14].
Recently, many physical, mathematical, financial, viscoelasticity, mechanical, and engineering problems have been described by the fractional variable-order pseudodifferential operators, see [24][25][26][27][28][29][30][31][32]. Also, the underlying differential operators in many dynamic processes appear as fractional and are dynamic in a sense that their derivative order is field-changing, i.e., it may change with space and/or time. Subsequently, a lot of scientists have been studying the characteristic of the variable-order fractional derivatives. In 1993, Samko with others [31] introduced this motivating expansion of the standard fractional derivatives by letting the order of the fractional derivative be not a constant but a function of an independent variable. After that, many authors have proposed various definitions of variable-order fractional derivative to fit coveted aims. Lorenzo et al. in [33] introduced in details the connotation of variable-order derivative and inspected some variable-order definitions. Coimbra [34] suggested a new definition for the variable-order fractional operator by using the Laplace transform of the Caputo fractional derivative.
The nonstandard finite difference methods (NSFDMs) were suggested by Mickens [35][36][37] to improve the discretization of some terms in the studied differential equation. Using a specific discretization and the denominator function, this approach will be more stable and more accurate than the standard approach [38,39], such that the proposed approach is not difficult to build [40]. The workable implementation of the NSFDMs exists in the research areas of chemistry, physics, and engineering [41][42][43]. In particular, the most fascinating usage is in mathematical ecology and mathematical biology [44,45], where the advantages of the NSFDMs have been shown markedly. Also, the active keeping property of the NSFDMs is well behaved in solving systems of fractional-order differential equations, like the fractional-order neuron model [46],the fractional-order Rössler model [47], and the Hodgkin-Huxley system with fractional order derivative [48].
In addition, the weighted average finite difference method (WAFDM) [49][50][51] can depend on the weight factor, and the implicit or explicit method (more stable or easy for coding, respectively).
In this work, we merge between those two techniques to find numerical approximating scheme for a time-space variable-order fractional nonlinear system of coupled hyperbolic and parabolic differential equations. This scheme has the advantages of these two approaches and therefore will has a larger stability region and will be more accurate. The introduced method is called the weighted average nonstandard finite difference method (WANSFDM) (see [18,50], and [52]).
Consider the system of general time-space variable-order fractional nonlinear, coupled, hyperbolic and parabolic partial differential equations in one-dimension with additional terms of a heat supply and a force term in the equation of heat conduction and in the equation of the motion, respectively: with initial conditions: and h 1 , are given smooth functions and the subscripts denote partial derivatives.
Many articles like [1,[53][54][55] have widely studied this system (1). When μ(x, t) = 2 and α(x, t) = 2, in [3,53] we found specific details of the meaning and application of such this system in physics. Also, in [50] the authors studied this system in case of time fractional derivative only. These works depend on the fundamental researches in which the theory of fractional thermoelasticity was introduced (see [56][57][58], and [59]).
In reality, the Caputo variable-order fractional derivative operator has many advantages when used to describe the derivative in initial value differential equations. The most important feature of variable-order Caputo's definition is that the conditions of the initial value for the variable-order fractional differential equations with the variable-order Caputo derivatives will be taken as in the case of differential equations of integer order, therefore, the time variable-order fractional derivatives are often taken in the Caputo sense. Also, it is usually to define the space variable-order fractional derivative as the Riesz variable-order fractional operator. The author in [13] concluded that "the complete theory of fractional differential equations, especially the theory of boundary value problems for fractional differential equations, can be developed only with the use of both left and right derivatives. " So, in this paper all the spatial derivatives are Riesz variable-order fractional derivatives, which include the right and left variable-order Riemann-Liouville fractional operators.
The fundamental aim of the present article is to carry out numerical treatment of the time-space variable-order fractional differential model (1) which arises in fractional variable-order thermoelasticity using the WANSFD technique. In reality, there are no studies in the literature using WANSFD technique to solve similar systems with variableorder fractional derivatives, so our attention to such an approach is motivated.
The article is arranged as follows: In the next section, we remind the elementary facts of the used methods (NSFDMs) and present many definitions of the derivatives of fractional variable order. In Sect. 3, in three different cases, we propose WANSFD schemes for the time-space variable order fractional model (1). We analyze, in Sect. 4, the stability of the suggested schemes and discuss their truncation error. Section 5 is devoted to some numerical consequences which are reported to exhibit the accuracy and the efficiency of the offered technique. Finally, in Sect. 6, we give conclusions of this work.

The nonstandard finite difference technique
The NSFDMs were suggested firstly by Mickens [35][36][37]. These techniques construct an approximating discrete schemes for partial differential equations (PDEs) or ordinary differential equations (ODEs). NSFDMs can preserve the properties of the exact solution of the studied PDEs or ODEs depending on the following principles [36]: 1. In general, the denominator functions for the discrete derivatives must be expressed in terms of the step sizes. 2. The orders of the discrete derivatives should be equal to the orders of the corresponding derivatives of the differential equations. 3. The nonlinear terms should be approximated in a nonlocal approximation manner. 4. The schemes must not have solutions which do not coincide with solutions of the studied differential models. 5. Specific conditions which keep the analytic solutions of the studied differential equations must also be specific discrete conditions for a scheme of the finite difference. In brief, for approximating dy dt using the Euler approach, we usually depend on y(t+h)-y(t) h but here we will use y(t+h)-y(t) is a continuous function of step size h which satisfies the following conditions: Furthermore, if there is a nonlinear term in the studied differential equation, we replace it in a nonlocal manner as, for example, x n+1 y n , x n y n+1 .

Fractional variable-order calculus definitions
In the literature, there are a lot of definitions of the operators of fractional variable-order derivatives (see, e.g., [60][61][62][63]). The space-fractional variable-order derivatives are usually given by the Riesz fractional variable-order derivatives. As for the time-fractional derivatives, they are usually defined in the Grünwald-Letnikov, Riemann-Liouville, or Caputo sense.

Definition 2.3 ([28])
The right-and left-sided variable-order Grünwald-Letnikov fractional derivatives are given by where [x] denotes the integer part of x and h is the step size.

Structure of the WANSFD for the proposed system
In the current part of this paper, we introduce WANSFDM for obtaining a discrete scheme for the system (1). Here, it is valuable to remind that the distinction between "the weighted average nonstandard finite difference technique" and "the weighted average standard finite difference technique" is analogous to the distinction between the nonstandard finite difference technique and the finite difference technique. The problem of finding the numerical solution of the system (1) is in approximating the fractional variable-order Caputo and Riemann-Liouville variable-order derivatives by WANSFD scheme when σ ∈ [0, 1] is the weight factor. Let N, M ∈ N and consider the coordinates of the mesh points, where We will denote, at a grid point (x n , t m ) = (nh, m t), the approximated values of u, θ by u m n , θ m n , respectively. The nonstandard difference approximation of variable-order Caputo derivative is introduced in the form: Similarly, Also, the approximations of variable-order Riemann-Liouville derivative using shifted Grünwald-Letnikov variable-order fractional derivative [29] are given as follows: where g α m n j = (-1) j α m n j .
In this paper we study the proposed model in three different cases: [1]) such that the constants have specific orders as follows: therefore system(1) can be written as follows: depending on the following difference approximations then system (10) takes the following form: The previous substitution gives an error, truncating error, symbolized here by T m n . Its estimation will be discussed in Sect. 4.2. After dropping this truncating error, we have the following computable difference scheme: and h 1 (x, t) = 0 (as in [55]) so model (1) has the following form: and then, using the difference approximations (11a)-(11f), system (14) can be written as follows: Neglecting the truncation error, then we have and h 1 (x, t) = 0 (as in [12]) so model (1) has the following form: where, depending on the difference approximations (11a)-(11f), model (17) then takes the form: All of the above schemes (13), (16), and (19) with their particular initial conditions and boundary conditions form an algebraic system of 2(N + 1)(M + 1) nonlinear equations with u m n , θ m n (n = 0, 1, 2, . . . , N , m = 0, 1, 2, . . . , M) as unknowns. The algebraic system can be solved using Newton's iteration method [64].

Analysis of the stability and truncation error 4.1 Stability analysis
To discuss the stability of the proposed schemes, we use the von Neumann-type stability technique and, because this technique deals only with linear problems, we have to linearize these schemes by supposing that the terms a 1 , b 1 , c 1 , d 1 , h 1 are constants and considering f = l = 0 (free source term) (see [1] and [65] for more details). To start our study analysis, we write system (1) in another form as follows: Letting u t = v, then we can express: Firstly, we write down system (21) in matrix form as follows: where Now, we write system (22) using WANSFDM as follows: which takes the following form: wherē w j = w j+1w j , j = 1, 2, . . . , m -2,w m-1 = -w m-1 , and w j = (j + 1) 2-μj 2-μ , j = 1, 2, . . . , m -1. if j > n + 1.
After some simplification we have: this scheme takes the following composite form: Depending on the von Neumann stability method, we presume such that β ∈ R, i = √ -1, Υ ∈ R 4×1 , and ξ ∈ R 4×4 is an amplification element. Putting Eq. (28) into Eq. (27), after some manipulation we obtain the following form of the amplification matrix: The scheme will be stable when ξ ≤ 1. As in [50,51,66], we consider on the right side the time-independent ξ = I so: Here ξ depends on μ m n , σ , h, and t, as we see from Eq. (30). Scheme (27) will be stable as soon as the following relation is satisfied:

Truncation error
Depending on definition of truncation error (15), (18), and (12), and from relations (11a)-(11f), we have Remark 4.1 It is well known that the Lax-Richtmyer equivalence theorem says: for consistent numerical approximations, stability and convergence are equivalent in the integer order differential equations. Yuste and Murillo in 2012 generalized this theorem to the fractional derivative case. They proved in [67,Sect. 4] that for a whole class of the fractional difference algorithms the consistency and stability of the difference scheme implies its convergence. Therefore our scheme is convergent under condition (31).

Numerical simulations
In f (x, t) = e -x cos t + (β 1 -1)e -x (1 -cos t) l(x, t) = (1 + a)e -x sin t -1 + be -x sin t e -x (1 -cos t) + 2α * e -2x (1 -cos t) 2 , and the exact solutions when μ(x, t) = α(x, t) = 2 (derivatives of integer order) are     Table 1 we see that E 1 and E 2 are reduced when N and M grow with σ = 0, 0.5 and 0 < t ≤ 1 and deduce that the approximations obtained by the explicit NSFDM do not converge because the condition of stability is not satisfied in this case.
Also, when 0 < t ≤ 5, we deduce the same conclusion from Table 2.
In Fig. 1 we see the approximation of u which was achieved by the introduced WANSFD scheme for various values of α, μ where σ = 0.5 and 0 < t ≤ 1.
In Fig. 2 we see the approximation of u which was achieved by the introduced WANSFD scheme for various values of α, μ where σ = 0.5 and 0 < t ≤ 5.  (14) for l(x, t) = -e -3t cos 3 x, and the analytic exact solutions for μ = 2, α = 2 (integer order derivative) are given by The boundary conditions are and the initial conditions are From Table 3 we see that E 1 and E 2 are reduced when N and M grow when σ = 0, 0.5, 0 < t ≤ 1 and deduce that the numerical approximation utilizing the explicit NSFDM is not convergent because the condition of stability is not satisfied in this case.
In Fig. 4 we see the approximation of u which was achieved by the introduced WANSFD scheme for various values of μ, α where σ = 0.5.
In Fig. 5 we see the approximation of θ which was achieved by the introduced WANSFD scheme for various values of μ, α where σ = 0.5. Figure 6 illustrates the behavior of the approximations of u obtained by the WANSFD technique when (σ = 0.5) at t = 0.5, for various values of μ and α = 2. Figure 7 illustrates the behavior of the numerical approximations of u obtained by the WANSFD technique when (σ = 0.5) at x = 0.7, for various values of α and μ = 2. Figure 8 illustrates the behavior of the numerical approximations of θ obtained by the WANSFD method when (σ = 0.5) at t = 0.5, for various values of μ and α = 2. Figure 9 illustrates the behavior of the numerical approximations of θ obtained by the WANSFD method when (σ = 0.5) at x = 0.7, for various values of α and μ = 2.
From Table 4 we see that E 1 and E 2 are reduced when N and M grow for σ = 0, 0.5, 0 < t ≤ 1 and deduce that the numerical approximation depending on the explicit NSFDM is not convergent because the condition of stability is not satisfied in this case.
In Fig. 10 we see the approximation of u which was achieved by the introduced WANSFD scheme for various values of α, μ and σ = 0.5.
In Fig. 11 we see the approximation of θ which was achieved by the introduced WANSFD scheme for various values of α, μ and σ = 0.5. Figure 12 illustrates the behavior of the numerical approximations of u obtained by the WANSFD technique when (σ = 0.5) at x = 0.3, for various values of α and μ = 2. Figure 13 illustrates the behavior of the numerical approximations of θ obtained by the WANSFD method when (σ = 0.5) at x = 0.3, for different values of α and μ = 2.

Summary and conclusions
In the paper, we introduced numerical simulations for nonlinear systems of time-space variable-order fractional, coupled hyperbolic and parabolic partial differential equations in one dimension utilizing WANSFDM. Analyzing stability of the proposed scheme is presented by a kind of John von Neumann method, also the relation which determines the truncation error is given. Some test examples are also presented. We compare the numerical results obtained by the used technique with the analytic exact solutions of the proposed model in the case of the standard derivative (integer order). We introduce many figures to show the behaviors of the solutions when μ and α are changed.